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Synopsis 

of the thesis entitled 

NUMERICAL TREATMENT OF SINGULARLY 
PERTURBED TWO-POINT BOUNDARY -VALUE SYSTEMS 

by 

Arindama Singh 
Department of Mathematics 

Indian Institute of Technology, Kanpur-208016, India 


Perhaps the human beings form a distinguished 
species of the genus 'Animals', because they possess 
the power of abstraction. It is this power of abstr- 
action that induces them to prefer mathematical 
models to describe real life problems. Then they 
analyse the models to understand and predict real 
life situations. 

Very often it happens that the models contain 
small parameters which one would otherwise like to 
neglect, (i.e., replacing small parameters by zero). 
It is all well if one is justified to do so. What 
if one is not justified in doing so? This is the 
beginning of singular perturbation phenomenon. 

Assume, for simplicity, that P_ is a problem 

O 

e-x pressed by ordinary differential equations, subject 



to appropriate boundary conditions, where only one 
small parameter e is involved. P is a singular 
perturbation problem, if the solution S 6 of P e 
does not converge (uniformly) to the solution S Q 
of the reduced problem P Q obtained from P g by 
letting e approach zero. However, in some cases, 
corrections to S Q could be done so that it approx- 
imates S g within, say, O(s^) accuracy. Asymptotic 
series solutions and matched asymptotic expansion 
techniques were found to be somewhat accurate in this 
direction. But the difficulties in using such tech- 
niques are well known; they lack simplicity and are 
very demanding on the part of the user. 

Numerical analysis accepted this new challenge. 
Consequently , numerical methods have been devised 
for special classes of singular perturbation problems. 
Beyond these special classes of problems most of these 
methods do not work. Methods which fulfil this need 
of transcending the subclasses of problems are neither 
simple nor cheap. 

In this thesis a general purpose method for 
solving systems of singularly perturbed two-point 
boundary-value problems in ordinary differential 
equations is proposed, analyzed and illustrated. 
Methods are also prescribed for some subclasses of 



this problem class , keeping the above requirements 
implicit. The thesis comprises seven chapters. 

Chapter 1 of this thesis presents motivation, 
a brief survey of the asymptotic and numerical 
analysis of such problems and summary of the present 
work. Results from various sources that are used 
in later chapters are also compiled to facilitate 
understanding the materials of this thesis. 

In Chapter 2, a simple partial decoupling 
transformation is constructed which decouples the 
slowly varying component from the original linear 
singularly perturbed system. Euler's Implicit 
(forward) and Explicit (backward) methods are shown 
to be accurate (of 0(h) or 0(h)) for the rapidly 
varying component. Non-commutativity of the differ- 
ence schemes and the decoupling transformations is 
also discussed. 

Based on the results of Chapter 2, an algorithm 
for singularly perturbed linear systems of two-point 
boundary-value problems is presented in Chapter 3. 
The effectiveness of the algorithm is demonstrated 
by some illustrative examples. 

tn Chapter 4, two special methods, namely. 
Boundary Value Technique and Cutting Point Technique 
are proposed. These techniques are applicable when 



the problems considered involve no turning points. 
These methods exploit the boundary-layer thickness, 
an indirect estimate of which is also derived in 
this chapter. 

In Chapter 5, the methods of preceding chapters 
are applied to non-linear singularly perturbed two- 
point boundary-value systems via a modified version 
of Newton's iteration. A model example of non-linear 
problem is solved to demonstrate the effectiveness 
of the methods. 

In Chapter 6, these methods are applied to 
optimal control problems both linear and non-linear. 

Chapter 7 provides a conclusion to this thesis 
and indicates scopes for further research. 

It is observed that the numerical methods 
presented have been found to be efficient over the 
conventional methods. The methods are at the same 
time, conceptually simple and need a very modest 
amount of problem preparation for computational 
purposes. Model problems of linear and non-linear 
types have been solved and numerical results are 
presented wherever necessary. The predicted 
accuracy could always be achieved without much 
computational effort. 



All the numerical results presented in this 
thesis have been computed on DLC-1090 Computer 
system available at the Indian Institute of Techno- 
logy, Kanpur, India. 



Chapter 1 


INTRODUCTION 


1.1 Introduction 

It is a fancy of pure mathematicians to invent 
and do such mathematics which has no obvious applica- 
tions: mathematics for the sake of mathematics. How- 
ever, most of mathematics that we do today have its 
origin in real life problems, or say, applications. 
Singular perturbation came into existence, not by 
fancy, but by necessity, thereby it belongs to the 
latter category. Ever since Prandtl’s [96] work in 
the beginning of this century, singular perturbation 
techniques have been a traditional tool of fluid dyna- 
mics. These techniques entered into various other 
areas of application, where ofcourse, the same ter- 
minology of ’boundary layer’, 'interior layer’, ’outer’ 
and ’inner’ expansions were already in use. Here, just 
to arouse one’s curiosity, we might mention its use in 
various applied areas such as fluid dynamics, plasticity, 
chemical reactor theory, nuclear reactor theory, plasma 
physics, aerodynamics, meteorology, rarefied gas dyna- 
mics, diffraction theory, reaction-diffusion process, 
non-equilibrium and radiating flows and etc. Later, 
we will consider an application area of singular pertur- 
bation : optimal control theory. Singular perturbation 



problems that cropped up in these and many other 
areas were treated almost exhaustively by asymptotic 
analysis. A breakthrough occured when numerical 
analysts turned their attention towards these problems. 
It was a challenge to numerical analysis. (We will 
shortly see why. ). Unlike the first expectation to 
see these two approaches as water-tight compartments, 
it turned out that wherever they are simultaneously 
applied, they rather stand complementary to each other. 
This combined effort turned out to be an efficient 
mathematical aparatus in investigating singular per- 
turbation phenomenon in various areas of engineering. 

The purpose of this chapter is to introduce and 
describe the problems those are treated in this thesis. 
Also a brief survey of the literature is presented. 

It should not only point out some of the problems that 
were or are being treated but also the manner in which 
these are being treated. Only a summary of some recent 
methods is presented. Obviously, a selection of tech- 
niques is done implicitly to present this summary, 
which is in some way relevant to the work pursued in 
this thesis. The methods so far developed, both pres- 
cribe as well as proscribe the approaches that should 
or should not be adopted, by recognizing the. problem 
classes. The problems at hand are such that, even an 
unskilled eye would see their worth as to be paid such 



a singular attention. We present here two examples 
to justify the singular attention and having a peep 
into the main difficulties involved. Without much ado, 
let us hear the monologues of first two explanatory 
characters of our drama. 

Example 1 . ([32]) 

Consider the following system of two-point boundary 
value problem : 

Y x = Y 1 + Y 2 (1.1a) 

0 < t < 1 

sy2 ~ ~~ y 2 ( 1 . lb ) 

with boundary conditions 

y x (o) = a , y 2 (i) = § (1.2) 

where ' . ' denotes d/dt, oc, j3 are some fixed real 
numbers, and e is a small positive parameter 
(0 < e << 1). 

Before going to any of the numerical methods , 
it is advisable to observe the equations carefully. 

This could best be done if we could get the exact 
solution. Fortunately, we do have the exact solution, 
i . e. , 

y-^tje) = (a + p eexp(l/e)/( 1+e) ) exp(t) 

- M exp((l-t)/e)/(l+e) 
y 2 (t,e) = p exp((l-t)/e) 


(1.3a) 

(1.3b) 



Now it is clear that changes very rapidly 
everywhere on (0,1] whereas y ^ blows up only in a 
small neighbourhood of t = 1 depending upon how 
small e is. This is a typical singular perturbation 
phenomenon, though the name 'Singular Perturbation* 
was given due to some other reason (cf. §1.2). How- 
ever, for numerical methods to be efficient, at least 
to solve such problems accurately, it is a warning 
that standard two-point boundary -value techniques 
with uniform grid should never be used, unless the 
grid length is very small (then of course, the round- 
off errors could surpass the accuracy level); and in 
all cases, the 'small neighbourhood* where y^ blows 
up should receive special care. Observing (l.l) 
through a quantitative eye, we find that the failure 
of the standard techniques can be attributed to the 
difference in order of magnitudes of the eigenvalues 
of the system matrix. The eigenvalues are 1 and -l/e ; 
when e<<l, this difference in order of magnitudes 
will lead to numerical instability of the standard 
methods. In fact, this is one of the main difficul- 
ties involved in such problems. Intentionally, to 
give a chill in one's spine, we have chosen such a 
simple example, where the system matrix is time- 
independent. In the next example, we take the system 
matrix to be time-dependent and enter into a bit dirty 
(but not unaesthetic) water. 



Example 2 . ([58]) 

Consider the following second order two-point 
boundary-value problem : 


ey + ty - y = 0 -1 < t < 1 (1.4) 

y(-l) - 1, y ( 1) = 2 (1.5) 

The exact solution of (1.4-1. 5) is : 

y(t,e ) = Cjt + c 2 (exp(-t 2 /(2e) ) 

t 0 

+ (t/e) f exp(-s 2 / (2e) )ds ) (1.6) 

-1 


where 

c 1 = -1 + 3 exp(-l/(2-e))/A 

(1.7) 


c 2 = 3/A 

(1.8) 


A =2 exp( -1/(2 e) ) + 1(e)/ s 

(1.9) 

with 

1 O 

1(e) = / exp(~s 2 /(2e))ds 

-1 

(1.10) 

Note 

that 1(e) = f2ne+ T.S.T. 

(1.11) 


where T.S.T. is an abbreviation for ' transcendent ally 
small terms'. With the help of (l.ll) the solution 
y(t,e) can be simplified as 


t < 0 
t « 0 


y 


-t - T.S.T. 
3fe/T2%J + T.S.T. 
-2t + T.S.T. 


t > 0 


( 1 . 12 ) 



This describes a ’corner’ at t = 0. Here again, 
we note that the eigenvalues of the system matrix are 
no longer time-independent as in Example 1. In fact, 
as t varies, the eigenvalues change their orders of 
magnitude on the interval of integration [-1,1], This 
is another main difficulty in the presence of which 
numerical stability of the standard techniques is 
jeopardised . 

The troublesome region in Example 1 is situated 
around the point t = 1, called a boundary layer. The 
situation in Example 2 is termed, in general, as a 
turning point (t = 0), which gives rise to interior 
or corner layers asymptotically. 

Naturally enough, the failure of standard tech- 
niques forces us to ask for better (specialized) 
techniques using non-uniform grids. We should not 
also forget to see from a user’s point of view. A 
user might not entertain himself to see the trouble- 
some regions beforehand. So we rather ask for a 
method or methods that indeed, succeed by finding out 
these troublesome regions simultaneously. This question 
will be answered in the following sections. 

There was a time when asymptotics had the monopoly 
in tackling such problems, and one had to seek an asym- 
ptotic solution unless one had already found out an 
exact (analytic) solution. So it is not advisable to 



neglect asymptotics altogether. Since our main 
emphasis is on numerical treatment, we will pass 
asymptotics over with a rapid pace. And in the 
later sections we will come across some of the 
numerical approaches that have been taken to tackle 
such difficult problems. 

1.2 Asymptotic Treatment of Singular Perturbation 

Let P be a problem in some differential model, 
which depends upon a positive small parameter e 
(0 < e <<l). Since e is small compared to 1, it is 
natural to consider also the degenerate or reduced 
problem P , obtained from P by replacing every 

O £ 

occurance of e in P £ with zero; here we assume that 
the problem P Q so obtained is meaningful. One then 
asks whether the solution S £ of P £ does converge uni- 
formly to the (reduced) solution S Q of the reduced 
problem P Q , as e-^O, If the answer is in the affir- 
mative, we call P , a regular perturbation problem, 

o 

whereas, if the answer is in the negative, we call P e 
a nonregular or singular perturbation problem. The 
word 'perturbation' is used viewing P Q as the original 
or unperturbed problem and thereby P , its perturbation. 

In the following we will be concerned with the 
study of the latter, namely, the singular perturbation. 

In particular we will restrict our attention to this 
singular perturbation phenomenon in ordinary differential 



equations and that to, only the two-point boundary- 
value problems (TPBVP). Unless and otherwise stated 
explicitly, singular perturbation problem will mean 
singularly perturbed two-point boundary-value problem. 

To be specific, we consider the problems of the 

type 

x = f(x, y, t , e) (2.1a) 

a < t < b. 

sy = g(x, y, t,e) (2.1b) 

with boundary conditions 

F( x ( a) , x(b) , y ( a) , y(b),e ) = 0 (2.2) 

where t£IR the real numbers, represents time, x(t,e) 
and y(t,e) are m and n dimensional real vector func- 
tions respectively; denotes d/dt? f and g are 

smooth real vector functions linear or non-linear of 
dimensions m and n respectively, and F in (2.2) repre- 
sents the boundary conditions so that (2.1) might be 
solved uniquely. In case, F is linear, it should 
represent (m+n) - linearly independent boundary 
conditions . 

In an excellent review article by O’Malley [85], 
history of asymptotic methods applied to such problems 
is pr-es-efrted. Starting from Prandtl's [96] paper on 
fluid dynamical boundary layers, which in fact, began 



the mathematical study of such problems, O'Malley 
lists many other important works on asymptotic 
methods. In asymptotics, a series solution of the 
form 

z(t,e) = z k (t).e , z = (2.3) 

is tried and under certain hypotheses such a series 

is found out for a specific member of the class of 

problems (2.1 - 2.2). The solution can be found out, 

. L, 

for practical purposes, upto 0(e) accuracy, for some 
finite L (typically L= 1 or 2). This rather concen- 
trates on the form of the solution, that too in a series 
form as in (2.3). However, since practically we 

9 

truncate the series in (2.3) after e or e , the terms 
left are represented as O(e^). This rather concentrates 
on the nature of the solution with respect to e . 

Kevorkian and Cole [58], for example emphasized 
on the nature of the solution on the whole interval 
[a,b], thereby studying the solutions from both the 
angles namely t and s . Sometimes, it happens that a 
series of the type 

CO . 

z(t,e) = 2 z k (t).(Ye) , z = 

k=o K 

approximates solution better than that in ( 2 . 3 ). Picto- 
rially, this represents the rapid change of the solution 
(see Example 1, § 1.1 . ) in a stretched independent 


(2.4) 



variable, which stretches t by a power of e. In 
(2.4) this power of e is taken to be Ve » the stretch- 
ing parameter. In a series solution given by (2.3), 
the stretching parameter is e itself. This idea can 
also be found in Eckhaus [27,28], though with more 
rigor. 

A simplified matching procedure has been pre- 
sented by Van Dyke [107] to tackle these problems 
efficiently. Now, in application, this matched asym- 
ptotic expansion process is used widely for such prob- 
lems. Also, a number of books, monographs and research 
reports have appeared which handle these problems mathe- 
matically. In a review article by Kokotovic [60], 
one can find an exhaustive bibliography on this topic. 
For an applied man new to this area, Ardema [4] pro- 
vides a good introduction. Before going to matching, 
it is worthwhile to look at Example 1 of Section 1.1 
once again. The second component of the solution 
vector, i.e., y 2 (t,e) blows up in a boundary layer 
at t = 1. Suppose we scale the equations (stretch 
the boundary layer) by a new independent variable 
a = (l-t)/s. Then y 2 (t, s) can be represented by 
y 2 (a, e) = p exp(cr) which does not blow up in a. Note 
that a scaling in equations produces also a scaling 
in the solutions. This region, the boundary layer 
region, is also referred to as the inner region; 



and the other region, beyond the boundary layers is 
termed as the outer region. In general, it is str- 
iking to take the stretched problem and solve it in 
the boundary layer region. This procedure is followed 
when the technique of matched asymptotic expansions is 
used. In the inner region we have the stretched prob- 
lem and in the outer region, the original one. But 
number of boundary conditions is still the same. To 
solve these problems, one extra boundary condition 
is imposed, which comes from the matching of two sol- 
utions of these two problems. The constants of inte- 
gration are determined by matching these solutions in 
an 'overlap region', i.e., where both the stretched 
and the unstretched solutions match. This matching 
is necessary since nothing about the exactness of 
the regions is known beforehand. These names such 
as inner region and outer region expansions, overlap 
regions etc., in fact, came from various application 
areas like fluid dynamics, optimal control, nuclear 
reactor theory etc. 

The results obtained, through matching coincide 
with those known by intuitive folkways of various 
applications. The most referred to texts include, 
for example, Bellman [13], Nay f eh [76,77], O'Malley [82], 
Eckhaus [27,, 28] , Vasileva [108], Tikhnov [105], 

Erdelyi [30], Ferguson [31]. 



The main problem with asymptotic techniques 
is that it is very demanding on the part of the 
user. The matching procedure is still not very 
simple to use. As Eckhaus [28] pointed out, one 
should know beforehand which stretching factor (some 
power of e , i.e., e , r: a rational number) is to 
be used. Our aim is to present numerical methods 
for users. For the fact that we do not miss a great 
name, this section is included. For more details 
regarding asymptotic expansion techniques, one 
should go through many of the literature not cited 
here. The key to those might be found out from 
O'Malley [85] and Kokotovic [60], 

1.3 Numerical Treatment of Singularly Perturbed 
Second Order TPBVPs. 

Leaving aside the analytic approaches or 
qualitative theories that have been applied to 
singularly perturbed two-point boundary value problems, 
we are forced to consider the asymptotics and the 
numerical means. Recently ,Non-~stand ard Analysis has 
been applied to singular perturbations [14,73]; but 
this only informs us about the geometrical behaviour 
of the solutions of a specific class of problems. 

For practical purposes, one needs quantitative means 
of solving the problems. In asymptotic analysis of 



such problems, the solution is represented as an 
additive sum of a certain (quantitative) function, 
and a term , the remainder or the residual term, 
which represents more of a property than a quantity 
[106, Ch-VIl], like O(e^) and etc. Numerical ana- 
lysis starts when we need to know the solution 
quantitatively and accurately. More precisely, 
while Numerical analysis provides us information 
regarding the solutions through quantitative - 
coloured glasses, Asymptotic analysis projects the 
solutions on a semi -quantitative , semi-qualitative 
tinted screen. However, numerical methods are in- 
tended at a broader class of problems and at the 
same time demands less from the user. Again, when 
for a specific problem one finds it very difficult 
to ride through asymptotics, it is natural to expect 
help from Numerical analysis. 

Unfortunately and not unexpectedly, the numer- 
ical analysis of singular perturbation is not that 
easy. As was remarked earlier, a non-uniform mesh 
is needed for the numerical schemes to produce accu- 
rate approximations. Note that this idea, if carried 
out might also reduce the number of grid-points, 
which is more or less necessary for an effective 
computation. Many such grid-generation processes 
have been used. Pearson [92], for example, used 



variable grid-length combined with a three-point 
difference scheme for the problems of the type 


e x + a(t)x + b(t)x = f(t) ~1 < t < 1 (3.1a) 

x ( — 1 ) = A , x ( 1 ) = B (3.1b) 


where denotes d/dt and & is a small positive 

parameter. Pearson's procedure in constructing a 
variable grid is conceptually simple. We require 
the grid to be dense in those regions where the 
solution changes rapidly. So the first step would 
be to know where does it really happen. For this 
purpose, he takes a uniform grid (with a fixed and 
modest e ) and solves the problem numerically with 
the help of a three-point difference-scheme. To 
test the rapid change in the solution, a tolerance 


limit 6 is prescribed for the difference 

df = - x jJ • d^ > § would mean refining the 

grid by introducing grid points between t. and 


ti+1* A smoothing process is also carried out to 
preclude abrupt changes within any mesh-interval. 

Then for smaller values of e , one takes the previous 


grid instead of a new uniform grid and then, repeat 
the whole procedure once again. Though this method 
is very simple to use, the major drawback is cost, 
i.e., in terms of machine-time. Pearson [93] also 
used this method for non-linear problems of the type 



F ( t , x , x , e x ) = 0 
x(0) = A, x ( 1 ) = B 


0 < t < 1 


(3.2a) 


(3.2b) 

The number of grid-points necessary for the numerical 
solution of these problems through this method is 
very large. For example, to solve a simple problem 

10~ 6 x + (x) 2 = 1 0 < t < 1 (3.3a) 

x ( 0) = x ( 1 ) = 1 (3.3b) 

upto five accurate digits, he had to construct a 
grid consisting of 4000 points. 

It may be remarked that a uniform grid is not 
expected to work satisfactorily due to the numerical 
instability of the problems considered. However, a 
fresh look at the numerical instability and an iner- 
tia towards using uniform mesh leads one to consider 
a directional difference scheme rather than the cen- 
tered difference schemes [24]. 

This idea was developed further by Dorr et 
al. [26], They applied a discrete maximum principle 
to obtain estimates on the numerical solutions of the 
equations (3.1). The choice of the direction in using 
the directional difference schemes (i.e., backward 
or forward) depends upon the sign of a(t) in (3.1). 
This also establishes the link with the asymptotic 
theory, as it is known that a boundary layer occurs 



at the left or right end point of the interval 
according as a(t) is negative or positive throughout 
the interval of integration [0,1], The main draw- 
back of this approach is that it can only be applied 
to a very restricted subclass of the classes of pro- 
blems (3. 1-3.2), particularly, where the maximum 
principle works. In the other direction, Dorr et 
al. [26] extended their results for the problem 
class which admits turning points and also to quasi- 
linear problems. 

Researchers sometimes feel reluctant to treat 
boundary-value problems and in particular, two-point 
boundary-value problems, because, in most cases, 
these can be transformed to initial value problems 
by the so called shooting methods. For example, 
Cohen [ 19 ] » Dorr and Farter [25] have considered the 
application of shooting methods for non-linear prob- 
lems of the type 

e x + f ( t , x, x ) x = 0 (3.4) 

and coupled equations 

x = f ( t , x, y) (3.5a) 

ey + g(t, x, x)y = c(t, x, x)y (3.5b) 


with appropriate boundary conditions. Kopell and 



Parter [64] analysed the problem 

ex = (x 2 - t 2 ) x (3.6a) 

x ( — 1 ) = A, x ( 1 ) = B (3.6b) 

again applying shooting methods. 

Numerical methods have also been devised 
which first approximate the co-efficients a(t), 
b(t) in (3.1) by polynomials and then apply some 
difference scheme. Ortiz [90], for example, esti- 
mated the errors in applying Tau method to the 
class of singular perturbation problems 

x + e"^p(t) x = 0 0 < t < 1 (3.7a) 

x ( 0) = 1, x ( 1 ) = 0 (3.7b) 

where 0 < e << 1 and p(t) is a polynomial approx- 
imation of a function not identically equal to 
zero on [0,l]. 

In asymptotics, the reduced problem plays 
an essential role. The solution of the reduced 
problem is called, as one expects, the reduced sol- 
ution. It happens that in the regions beyond the 
layers, the reduced solution represents or approxi- 
mates the original solution fairly well. Flaherty 
and O'Malley [35] have given a numerical algorithm 
for singularly perturbed second order TPBVPs , or 



more generally, for second order stiff TPBVPs, based 
on the reduced problem. The characteristic roots 
of the reduced problem enables them to obtain the 
boundary layer correction terms numerically. These 
correction terms are corrections to the reduced 
solution(s) in respective regions. The reduced 
solution could also be obtained from the reduced 
problem numerically. This is a beautiful example 
of the cross-breeding of the two approaches : 
asymptotic and numerical. 

Carried on by the asymptotic experiences , 

Lorenz [72] has considered the numerical solution 
of the problem 

ex = a(t,x) x + J3(t,x) = 0 0 < t < 1 (3.8a) 

x ( 0) = A, x ( 1 ) = B (3.8b) 

He splits the interval of integration [0,l] into two 
parts and then on one part solves the reduced problem, 
whereas on the other part, solves the original problem. 
The latter part being stretched appropriately before 
the difference schemes are applied. 

Recently, Kadalbajoo and Reddy [47-50] have 
given some methods for singularly perturbed two-point 
boundary-value problems of second order, which go in 
the same vein of observing the problem from both the 
asymptotic as well as numerical points of view. 



Also many non-asymptotic methods have been 
invented. For example, Roberts [97] has given a 
method, there called, Boundary Value Technique for 
singularly perturbed TPBVP's. He has also modified 
the method a bit so that it could be applied to 
problems of an extended class [98-99] . These methods 
have been found successful even in many non-linear 
cases . 

Among others, the three-point difference 
schemes of Kellogg et al . [57], the chopping proce- 
dures of Sakai [lOO] , Sakai et al. [lOl], exponen- 
tial tri-diagonal difference scheme of Berger et 
al. [15] , collocation methods using cubic polyno- 
mials and cubic splines by Flaherty and Mathon [34], 
box and trapezoidal schemes studied by Weiss [113] 
note-worthy. Axelsson and Carey [12] constructed a 
regular modified problem to discuss the boundary 
layers separately and derived thereby the estimates 
for boundary layers which are uniform in terms of 
the small (but immensely notorious) parameter e . 

1.4 Numerical Treatment of Systems of Singularly 
Perturbed TPBVPs. 

Compared to numerical treatment for handling 
syst-ems of singularly perturbed two-point boundary 



value problems much has been done by the asymptotics. 
In various treatments of such problems, several 
problem classes have been considered. Among these 
problem classes, much attention has been paid to 
linear homogeneous systems because of their occur- 
ence in a wide range of applications. For example, 
in Optimal Control Theory, when a scalar (quadratic) 
cost is optimized subject to the states modelled by 
ordinary differential equations containing a small 
parameter £ (0 < £<< l), which typically represents 

small mass, inductance, capacitance etc., it happens 
that the resulting differential equations to be sol- 
ved become singularly perturbed. Kokotovic [60] 
claims that almost all control problems are singu- 
larly perturbed; one does not recognize them because 
it is customary to neglect the small parameters. 
Though several such parameters might appear in a 
specific model giving rise to singularly perturbed 
systems with several parameters, we will be inter- 
ested here in describing only one-parameter cases. 

The linear homogeneous system of singularly 
perturbed TPBVP's can be represented as follows : 

x = Ax + By (4.1a) 

0 < t < 1 

ey = Cx + Dy (4.1b) 



(4.1c) 



where x(t,e), y(t,e) are vector functions of dim- 
ensions m,n • A(t,e), B(t,e), C(t,e), D(t,e) are 
given matrix functions of dimensions mxm, mxn, nxm, 
nxn; and (4.1c) represents (m+n) - linearly indepe- 
ndent boundary conditions? ' denotes d/dt , and e 
is a small positive parameter, 0 < e << 1. 

Asymptotic analysis of this problem (4.1) reve- 
als the fact that if, A, B, C, D are smooth bounded 
matrices, then x(t,e) is a slowly varying component 
of the solution and y(t,e) is the rapidly (fast) 
varying component. Under certain hypotheses the 
asymptotic solution has been obtained. To name a 
few, credits might go to Vasileva [108], Tikhonov [105], 
(see also [ill]) , O'Malley [82], Kokotovic et al. [62], 
Hopensteadt [44], Sannuti [103], Finden [33] ,Nayf eh[76] , 
Nipp [79] and for their other works in the related 
fields as well as that of their co-workers. However, 
for a detailed bibliography one can see Kokotovic [60] . 

Various assumptions have been taken in order to 
ensure that the solution z 6 (t) of (4,1) approaches 
the reduced solution 2 (t) of the corresponding 
reduced system ( z = j~ * J ) 

x = Ax + By (-4.2a) 

0 = Cx + Dy • (4.2b) 

with appropriate choice of m-boundary conditions from 



(4.1c), in the limit as e — > 0+ . Recently, 

Campbell [17] has given a new approach for initial 
value problems. Also, differential inequalities 
have been applied to guarantee the existence/unique- 
ness as well as this limiting behaviour of the 
solution z (t). 

, S 

Howes [45] considers singularly perturbed 
second order non-linear systems of the type 

ex+F(x)x + g(t,x)=0 a<t<b (4.3) 

where F is a matrix valued function. Sufficient 
conditions are given in [4b] in order that a solu- 
tion of (4.3) displays a layer behaviour, i.e., non- 
uniform dependence of the solution on t and e 
as e — ■> 0+ . 

Flaherty and O’Malley [35] combined the asymp- 
totic analysis and numerical analysis, in the sense 
that, to construct a numerical solution, they used 
the asymptotic information. Though, they have gen- 
eralized their method to quasi-linear scalar case, 
they have not extended it to systems. However, as 
we have remarked earlier, the asymptotic analysis of 
singularly perturbed systems of two-point boundary- 
value problems demand, also, a variable mesh const- 
ruction for the effectiveness of the numerical 
algorithms. This means, in different regions of the 
interval of integration, different grid-lengths 



must be used as well as a stretching of the small 
regions helps in achieving a comparatively flaw- 
less computation. Abrahamson et al. [l],for example, 
analysed the second order system 

ex + A(t) x + B ( t ) x = F(t) 0 < t < 1 (4.4a) 

x ( 0) = a, x(l) = p. (4.4b) 

They assume that 

A(t) = diag( A 1 ( t ) , A n (t)) (4.5a) 

A = A* (4.5b) 

where A* is the adjoint of A, and there is some con- 
stant r) > 0 such that 

A 1 < r] < 0 and A 11 > r\ > O (4.5c) 

for every t, C < t < 1? where e.g. < r\ means 

T 

that all the entries in A are less than or 
equal to r). 

Introduce the notation 



r- — 


r~ _ ~t 


x 1 


F 1 

X = 

x 11 
— • — * 

, F = 

pH 


corresponding to the partition of A. This partition 
plays its role in letting the boundary layers appear 
in narrow regions at x = 0 and at x = 1 for the 



II I 

variables x and x respectively. This has 
been proved by finding an a priori estimate of the 
solution of (4.4). In constructing numerical algo- 
rithms it is important that the numerical algorithm 
should be able to work efficiently when the mesh 
size h >> e. Otherwise, one has to take h very 
small to carry out the computation thereby causing 
roundoff errors and resulting in prohibitive computer 
time . 

Abrahamson et al. [ 1 j used the following diff- 
erence approximation. The problem (4.4) is approxi- 
mated by 

(s + h C k )D + D_ u k + A k D o u k 

+ L 1r (B,u) = L 2 k (F) (4.6a) 

u 0 = « , U N = £ (4 * 6b) 

where C(t) = diag(C I (t), C I3: (t)) is a positive 
definite matrix such that C(t) >. -j|diag( -A F ( t ) ,A FF ( t ) ) , 
D , D and D are the usual forward, backward and 
central difference operators and, are linear 

operators of the form 

L lk (B - u) = B -lk u k-l + B ok u k + B lk u k+l 

L 2k (F) = C -1 F(t k - 2 } + C o F( V + C l F(t k + h/2) ‘ 



The coefficients B . . , C. ' s are chosen to make (4.6) 

**• J _L 

as accurate as possible; it would also make the 
boundary layers as sharp as possible. They have 
used Richardson extrapolation in the interior of 
the interval to achieve the accuracy upto 0(e). 

The smoothness of the numerical solution has also 
been shown except in the boundary layer regions of 
thickness of 0( e| In e |) and the accuracy achieved 
is 0(e), and the rate of convergence is 0(h). 

Following this approach as well as Pearson's [92- 
93] approach of variable mesh construction, Kreiss 
and Kreiss [67] developed a method for this class 
of problems. A variable mesh is suggested instead 
of fixing the mesh beforehand. Starting with some 
reference mesh, a difference scheme of the type (4.6) 
is used to compute a numerical solution. Then, a 
new mesh is constructed by adding or deleting mesh 
points according to the variation of the numerical 
solution. This novel feature enables them to extend 
the method to the cases when the condition A = A* 
in (4.5b) is not satisfied; as well as this method 
is also applicable to the non-homogeneous form cor- 
responding to (4.1). The assumptions (4.5) are con- 
nected with the two main difficulties associated with 
systems like (4.1) and its non-homogeneous form. The 



difficulties encountered are (see the Examples 1, 2 
in $1.1 ) : first, the eigenvalues of the system 

matrix 


A B 



vary widely in magnitude, and secondly, these eigen- 
values change their orders of magnitude when t varies 
along [0,l] (see [66, 75]). 


Methods have also been developed for a special 
class of problems, namely, singular-singularly per- 
turbed two-point boundary value systems. For example, 
O’Malley [84] considers systems of the type 


x = A(t,e)x + B(t,e)y + C(t,e) (4.7a) 

ey = E(t,e) x + F(t,e)y + G(t,e) (4.7b) 


(4.7) is called singular by assuming that 
the matrix F(t,0) is singular. A full set of 
asymptotic solutions is constructed assuming further 
that F(t,0) is block-diagonalizable and the reduced 
problem is consistent and a stability condition 
(defined in [84]) holds. This would give rise to 
numerical methods mentioned earlier, for example, as 
in [35]. 



Ascher [7] considers the following system of 
TPBVP's : 

ex = A(t,e)x + f(t,e) 0 < t < 1 (4.8a) 

B Q (e) x(0, e) + ( e) x(l, e) = j3( e) (4.8b) 

where x(t, e) is a vector valued function, and 
f(t, e) is a vector valued linear or non-linear 
function, A(t,e), a matrix valued function of appro- 
priate dimensions; and (4.8b) represents boundary 
conditions. The reduced problem corresponding to 
(4.8) is assumed to be singular. For numerical 
approximation, families of symmetric difference 
schemes are used. Convergence results for the 
regular singularly perturbed case are extended with 
reference to a modification in mesh selection. Con- 
vergence of mid-point rule and trapezoidal rule are 
proved. Ascher proved this convergence by assuming 
some conditions on the difference schemes rather 
than on the problem (4.8) itself. 

Recently Kreiss et al. [68] have given a method 
for stiff two-point boundary-value problems. This 
method uses a grid-generation process closely related 
to that of Abrahamson et al. [ 1 ] . The difference 
scheme is a combination of Euler's Implicit, Euler's 



Explicit and Trapezoidal rules. However, the grid- 
generation process consumes too much of machine 
time and it also involves a large number of function- 
evaluations . 

The non-linear singularly perturbed two-point 
boundary -value systems are treated with shooting 
methods, as is the case for non-linear stiff counter- 
parts. Using Newton's Iteration [89], the non-linear 
problems could be realized as the limit of a sequence 
of linear problems and then the methods for linear 
problems could be used [68]. 

Due to the undoubted status of shooting, mathe- 
maticians are reluctant to accept the adjective 'stiff' 
being applied to two-point boundary-value problems 
( 8.2 in [2]). However, the case of singular 

perturbation is non-controversial. 

Among others, difference schemes of Dorr [24-25], 
spline collocation codes of Ascher et al. [6-7, ll] 
could very well be incorporated for systems of general 
linear or non-linear singularly perturbed TPBVP's. In 
fact, the collocation codes were already constructed 
for this problem class, see e.g., Ascher and Weiss [8™ 
10], Centred difference schemes were analyzed by 
Kreiss [65] when applied to singular perturbations. 



An inclusion relationship between perturbed colloca- 
tion and ftunge-Kutta methods was investigated by 
Nj^rsett and Wanner [80]. Uniform stability of dis- 
crete and continuous singularly perturbed boundary- 
value problems was studied by Niederdrenk and 
Yserentant [78], Weiss [113] analyzed the box scheme 
and the trapezoidal rule for linear singularly per- 
turbed boundary-value problems. Smith [104] consi- 
dered a singularly perturbed boundary-value problem 
arising in the physical theory of semiconductors. 

For difficult topics like singular-singularly perturbed 
TPBVP's, Vasileva and Butuzov [110] have given the 
asymptotic solutions. The one-sided difference schemes 
'are shown to be applicable in nonlinear singularly 
perturbed problems by Osher [91]. 

1.5 Summary of the Thesis 

In the previous sections, we have come across 
various classes of problems as well as asymptotic and 
numerical methods for their solutions. Both the asym- 
ptotic and numerical methods aim at finding accurate 
solutions to concerned problem (classes). The accu- 
racy in asymptotic expansions is achieved upto an 
order in some power of e (0(e n )) and in numerical 
techniques, the accuracy is achieved upto an 



order in some power of the maximum grid length h 
(0(h m )). Next to accuracy, we also need a numerical 
method to be efficient with respect to cost, i.e., 
less consumption of machine time as well as less 
interference of the user. The other direction is 
to find out general purpose algorithms which would 
handle various classes of problems e.g., singular- 
singularly perturbed problems, singular perturbation 
problems with turning points and non-linear singular 
perturbation problems. 

Asymptotic expansion techniques, infact, as the 
literature shows, cover almost all these problem 
classes. However, for a naive user's point of view, 
this is a difficult task. For every particular pro- 
blem, he has to find out the inner expansion and the 
outer expansion and match them to determine the co- 
efficients of the expansion terms. This matching 
requires skill, insight into the concerned problem 
and a great deal of experimentation. An engineer or 
an applied mathematician needs efficient techniques 
which he can apply routinely in order to avoid deal- 
ing with methods which demand so much of skill, in- 
sight and experimentations. 

In response to this need for a fresh approach 

to such problems, some efficient and simple methods 
are proposed and illustrated in this thesis. 



In Chapter 2, a simple decoupling transforma- 
tion is constructed, and analyzed. Unlike the usual 
Lyapunov type transformation, this transformation 
partially decouples a singularly perturbed system. 

The s lowly -varying component is decoupled completely, 
whereas the rapidly varying component involves the 
slowly varying one. Euler's forward and backward 
implicit methods are shown to be accurate (in com- 
plementary cases) for the rapidly varying component. 
The non-commutativity of the decouplings and the 
difference schemes is also shown. 

In Chapter 3, an algorithm based on the results 
of Chapter 2 for linear systems of singularly per- 
turbed two-point boundary-value problems is presented. 
Numerical experimentations are carried out to demon- 
strate the effectiveness of the method and thus, the 
algorithm. 

Based on some crucial observations met in the 
process of demonstrating the effectiveness of the 
hybrid method, some other numerical methods are pre- 
sented in Chapter 4. These methods, namely, the 
Boundary Value Technique and the Cutting Point Tech- 
nique also take into account the facts known from 
the asymptotic analysis of problems, where the turn- 
ing points are absent. These methods also exploit 



the numerical boundary layer thickness, an estimate 
of which is also found out in this chapter. The 
effectiveness of these special methods are illustrated 
by examples. 

In Chapter 5 the methods are applied to non- 
linear problems via a modification of Newton's 
Iteration. A model example of non-linear singular 
perturbation is taken from Kevorkian and Cole [58] 
for numerical illustrations. 

In Chapter 6, these methods are applied to 
some optimal control problems, both linear and non- 
linear. 

Finally, in Chapter 7, a veil is drawn as to 
epitomise this work and scope for further research 
is also indicated. 

In a nut-shell, the numerical methods presented 
in this thesis for solving systems of singularly per- 
turbed two-point boundary-value problems in ordinary 
differential equations have been shown to be efficient 
over the conventional methods. The methods are con- 
ceptually simple and need a very modest amount of 
problem preparation for adapting these to computer 
implementation. Various kinds of linear and non- 
linear examples have been solved and numerical results 



are presented wherever necessary. It is observed 
that the accuracy predicted can always be achieved 
with practically no computational effort. All the 
numerical results presented in this thesis have 
been computed on DEC-1090 computer system at III 
Kanpur . 


1.6 Preliminary Results 

In this section we just introduce some results 
from various sources, that will be needed in the 
sequel. The theorems below will appear disconnected 
from each other; and no attempt is made to connect 
them. The purpose of adding this section is that it 
will facilitate presentation of the contents 
of this thesis. The results given below are non- 
trivial and they have no obvious connection to 
Singular Perturbation. Thus, when we refer to these 
results, it will be convenient for one who does not 
want to go through the original sources but realise 
the importance for their applications. 

The first theorem tells us about the block-diagonal- 
ization of a general matrix. 



Theorem 6.1 ( [ 38 ] ) 


Let I be a closed sub-interval of [0,l], Let 
the eigenvalues of a time dependent square matrix 
A(t) of order m be X^(t), i = l,2,...,m. Let 
<p(t,X) be the polynomial with spectrum 
{ X i (t) | i = 1,2,. . . ,m }. Let cp^ ( t ,X ), <P 2 (t,X ) 
be relatively prime polynomials with spectra 
{X i (t) i i = 1 , 2 ,..., r} and {X i (t) | i = r+l,...,m } 
respectively. Suppose that <p(t,^) = <P^( t , X) 4 <P 2 ( t , X) 
for all t £ I. Then there is a smooth invertible 
matrix M(t) of order m, t £ I, such that 
M“ 1 (t) A(t) M(t) = diag(A 1 (t), A 2 (t)) where the 
eigenvalues of A^t) and A 2 (t) are 

X x (t) ,. . . , X (t) , and * r+ 1 (t),..., *- m (t) respectively. 
Theorem 6.2 ( [ 68 ] ) 

Let 6 > 0, > 0 be constants and consider 

a grid function satisfying 
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The theorem below is the generalized Taylor’s theorem 
for operators. (Linz [70]). 


T heorem 6.3 


Let P : X — > Y be an operator between two 
Banach spaces X and Y such that P is n-times 
continuously differentiable in a neighbourhood 
N(x o ,r), r > 0 of the point x q in X. Then, for 
all x in the interior of N(x ,r). 


|P(X) - P(* 0 ) m , 

m=l 


i y £ S L(x ,x) II pW <V> H 
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where L (x Q ,x) is the line segment between x Q and 
x, and ||.|j denotes the supremum norm. 

Note that the expression P^ m ^ (x Q ) (x-x Q ) m is 
meaningful since P^ m -*(x 0 ) is an m-linear operator. 



P^^(x 0 )(x-x )^, for example, means 

((P^ (x Q ) (x~x 0 )) (x-x q )) (x~x Q ). 

The next theorem is the Pontryagin's Maximum 
Principle [95] whic-h can be applied to reduce 
optimal control problems to a system of differential 
equations . 

Theorem 6 .4 

Let the scalar cost be defined by 

h 

J = / x (x( t ) , u( t ) )dt which is to be maximized by 
t 

o 

selecting a control vector u(t) subject to the 
states x(t) satisfying the differential system 
dx/dt = f(x,u), x(t Q ) = x . The vector functions 
f(x,u) and the partial derivatives of each of its 
components are assumed to be given and are continuous 
on the direct product X x U, where U is the clo- 
sure of U in lR r , the control space and X the 

state space. Let x° satisfy the differential 
, o 

equation = f°(x( t).,u( t ) ) and u be a given 

line parallel to the x° axis which passes through 

the point (0,x(t^)). Let u(t) , t < t < , 

be an admissible control such that the corresponding 

trajectory x(t) which begins at the given point 

x at the time t is defined on the entire interval 
o o 



t <. t < , through a point on the given line %. 

In order that u(t) and x(t) be optimal it is 
necessary that there exists a non-zero absolutely 
continuous vector function ¥ (t) = ( ¥ (t),..., ¥ n (t)) 

corresponding to the functions u(t) and x(t) such 
that 

(i) the function H(¥(t-), x(t), u) 

n 

= I ¥ • f^(x,u) of the variable u in U 
3=o J 

attains its maximum at the point u = u(t) almost 
everywhere in the interval t < t < tj , 

(ii) at the terminal time t, the relations 

%(t x ) < 0, S u U ^ u ( ¥ (t x ), x(t 1 )) = 0 (6.1) 

are satisfied. 

Furthermore, it turns out that if ¥ (t), x(t) and 

o 

u(t) satisfy the system x = H^, , ¥ = - H x and 

condition (i) ; then the time function ¥ Q (t) and 
sup H(¥(t), x(t), u) are constant. Thus, the above 

u eu 

condition (6.1) might be verified at any time t, 
t <. t < tj_ and not just at t-^ . 

The next result, namely, the implicit function 
theorem is taken from Lang [69]. 



Let U, V be open sets in Banach spaces E,F 
respectively, and let f : U x V -—> G be a p-times 
continuously differentiable function, G a Banach 
space. Let (a,b) £ U x V and assume that 
D 2 f(a,b) i F *— > G is invertibel. (D 2 denotes the 
partial derivative with respect to the second argu- 
ment). Suppose that f(a,b) = 0. Then there exists 
a continuous function g : U ™~-> V defined on an 
open neighbourhood U q of the point a such that 
g(a) = b and f(x, g(x))=0 for all x in U Q . 

If U o is taken to be sufficiently small ball, then 
g is uniquely determined and also is p-times continu- 
ously differentiable. 

Next result needs the following definitions [89]. 

Let L(R n ,R rn ) be the vector space of linear operators 
from R n to R m (R^ being j-dimensional real space). 

Definition 1 

A mapping F : D C R n — > R m is G~dif ferentiabl 

at an interior point x of D if there exists a linear 

operator A£L(R n ,R m ) such that for any h£ R n , 

lim | | F(x+h) - Fx - t Ah | j =0. The G-derivative 
t-o- 

of F is denoted by F'. 



Definition 2 


The mapping F : D C R n — — -> R m is Frechet 
differentiable at x £ Int(D) if there is an 
A £ L (R n ,R m ) such that 

lim (1/j | h i 1 ) | |F(x+h) - Fx ~ Ah | | = 0. 

h-*o 

The linear operator A is called the Frechet- 
derivative of F at x. 


Definition 3 

Let I : x k+1 = G x k , k = 0,1,... .be an iter- 
ative method (process), where G : D C R n R n . 

Then x* is a point of attraction of the iteration 

I if there is an open neighbourhood S of x such 

o k 

that S C D and for any x £S, the iterates Cx } 

-&* 

defined above all lie in D and converge to x . 


Definition 4 

1c n 

Let {x }££ R be any convergent sequence with 
limit x* . Then the quantities : 




) = 


0, if x k = x’ for all 

but finitely 
many k. 

lim sup — 7 * -JJ , if x k ^ x* for all 

k~*°° j | x K — x^ | j p but finitely 

many k. 




CO 


otherwise 



Definition 5 


Let C(l,x*) be the set of all sequences with 
limit x* generated by an iterative process I. Then 


Qp(I,x*) = sup { Q p x* } I { x k }£ C ( 1 , x* ) }, 1 < P < «» 


are the quotient convergence factors of I at x ,r with 

lr 

respect to the norm in which the Q p { x } are computed 


Oq(I> x *) = 


co, if q (i, x *) = o, Vpe [i,~) 

inf { p £ [l,«>) j Qp( I »x*) = °° ^otherwise 


is the Q~order of I at x\ 


Definition 6 

lr 

Let { x JC R be any sequence that converges 
to x*. Then the numbers 


Rptxh = 


lim sup 

| x k -x* 1 1 1/k , 

if p = 1 

k-^co 



lim sup 1 

iik *ii l/k p 


! 1 x -x | | ' 

if p > 1 

k-» oo 


are the root convergence factors of the sequence 

{x }. If I is an iterative process with limit 
* * 

point x and C(I,x ) is the set of all sequences 
generated by I which converge to x , then 



are the R -factors of I at x*. 
The quantity 


0 R (I, x*) = 


03 , if R p( I ?xi ' > = °» V P£ 

inf { p£ [1,~ ) |R (I,x*) = 1 } 5 otherwise •, 


is called the R-order of I at x*'. 

Now we state the convergence theorem of the Newton's 
iteration [89, Thm 10.22]. 


Theorem 6.6 

Assume that F : D C R n — R n is G-differen~ 
tiable on an open neighbourhood S Q <Z D of a 
point x* £ D for which Fx* = O, and that F' is conti- 
nuous at x* and F'(x*) non-singular. Then x* 
is a point of attraction of the Newton iteration 

I : x k+1 = x k - F'(x k )~^ F x k , k = 0,1,... and 
R^(I,x*) = Q X (I, x*) = 0. If, in addition, there are 
constants a < °o and p£(0,l] such that 
| |F ' (x) - F' (x*) | | < a |x-x*| p , V x £ S Q , then 
0 R (I,x*) _> Oq(I,x*) _> 1+p. Finally, if F is conti- 
nuously Frechet-diff erentiable on S Q and the second 
Frechet-derivative of F exists at x* and satisfies 
F" ( x* )h 4 0, Vh £R n , h 4 0 then 0 R (I ,x*)=0 Q (I,x*)=2 ; 
where the last F" denotes the second Frechet deri- 


vative of F. 



Chapter 2 


A METHOD FOR LINEAR SINGULARLY PERTURBED SYSTEMS 

2.1 Introduction 

For numerical integration of ordinary differen- 
tial equations , the first and foremost task is to 
find out an appropriate grid. This appropriateness 
can be best represented via the assumptions on the 
coefficient matrix and the non-homogeneous terra, if 
any, occuring in the differential equation. This is 
exactly done, for example, in Kreiss et al. [68] who 
prove the boundedness of the solution and its deriva- 
tives upon assuming some numerically verifiable condi- 
tions on the system matrix and the non-homogeneous 
term for stiff linear systems. Our concern here is 
a special subclass of these stiff systems, namely, 
the singularly perturbed systems. Ferguson [32] 
considers the two-point boundary-value systems of 
the type 


"’® Yl~ 

O 


A 11 A 12 

6 A 13 


y l 


1 f l- 

y 2 

• 

= 

A 21 

A 22 

A 23 


y 2 

+ 

f 2 

e y 3 


& A 0 -j 

A 32 

A 33 - 




„ f 3„ 


(1.1a) 

with the b-oundary conditions 
B q y(O) + B 1 y( 1) = 3 2 


(1.1b) 



where y 


*1 

y 2 

y 3 


, ( 1. lfc?)represents n-linearly independent 


boundary conditions, y £ IR n * 0 < e << 1, a para- 

meter, ' . ' denotes d/dt. Assuming that the matrices 
A. .(t), f. (t) are smooth (possess continuous deriva- 

J J- 

tives upto some finite order) functions of t satis- 
fying the stability conditions 


Re X(A U ) < < 0 < n < Re\(A 33 ) 

he proves the uniqueness of a solution y which 
admits of the following a priori bounds : 


llyxll 

< K(| Yl (0)| + jy 2 (0)j + e ( y 3 ( 1 ) | 

+ s llfilt) 

i=l 1 

(1.2a) 

lly 2 ll 

< K(ejy 1 (0)| + |y 2 (0)| +ejy 3 (l) 

1 * ^ llfill) 

i=l 1 

(1.2b) 

lly 3 ll 

< K(e| Vl (0)| + ! y 2 (0) I +e|y 3 (l) 

1 + 2 W^W 

i=l 1 

(1.2c) 

with 

1 

| |g| I =f |g(t)jdt, and K(> 0) 

is a real 


o 


constant. The results corresponding to (1.2) for the 
general class of linear singularly perturbed two-point 
boundary value systems : 


O' 

- y r 


A 11 A 12 


" y- 

+ 

f i 

e°V 2 


; A 21 A 22 


y 2 


f 2 


, 0 < t < 1 

(1.3) 



with similar boundary conditions (1.1b) are also 
given in his thesis [3l], 

In addition to these a priori estimates of the 
solution, we also know the geometric behaviour of 
the solution from the asymptotic analysis. In 
general, the boundedness results are proved under 
the conditional stability criterion, that the eigen- 
values of the system matrix A(t£ ) = 

are bounded away from the imaginary axis. Under this 
conditional stability criterian, the rapidly varying 
component y 2 (t,e) of the solution y(t,e) can 
become a combination of two pieces contained in dif- 
ferent manifolds. The dichotomies (ordinary or 
exponential) present in the fast varying component 
gives the motivation to decouple the system first 
into two modes: slowly varying and the rapidly varying 
one. From Coppel [20-21] the theorem below follows $ 
(Kokotovic et al. [62, Ch-5] ) « 

Theorem 1,1 

Let B(t) be a continuously differentiable . 
matrix function such that for every t£ [0,l], 

I I B ( t ) | | < M, | | a ( t ) j | < 6 , 


A 11 A 12 
A 21 A 22 


and 



there exists a > 0 : Re X(B)(t)) < - a . (1.4) 

Then there exists e 0 > 0 such that, for 0 < e <e Q 
the system e z = 3(t)z has a fundamental matrix 
X(t) satisfying the inequality 

I|x(t) X-tsJll < K e- a(t - s) (1.5) 

for t >_ s , for some const. K > 0. Moreover, if 
condition (1.4) is replaced by 

there exists a > 0 : Re \(B) >_ cc , (1.6) 

then the inequality 

| |x(t)X~ 1 (s)| | < K e ~ a(s ~ t ) (1.7) 

for s >_ t, holds for some constant K > O. 

Remark 1.1 

The constants e Q and K, in case (1.6) is 
satisfied (and thus in (1,7)) might differ from those 
for the case where (1.4) is satisfied (and accordingly 
in (1.5)). For the second case (i.e., (1.6)) we write 
the corresponding and K. by and K Q 

respectively. 

In the following we abbreviate the above matrix 
product X(t)X ml (s) by refering it to as 'the matrizant' 



<p(t,s) = X(t) X _1 (s) ; and 

<p(s y t) = X(s) X" 1 (t). As such the inequalities 
(1.5) and (1.7) will be interpreted as giving 
estimates of the matrizants. 


2.2 A Decoupling Transformation 

The usual way of decoupling the system 


*1 


A ll (t) A 12 (t) 


— 


h(t)' 


rr 






f y 2j 


A 21 (t) A ?2 (t) 

! 

_ V 2_ 

i 

_ f 2^ t ^ _ 


for 0 < t < 1, is done by finding a decoupling 

r I eH(t) 

~ m 


Lyapunov transformation T(t) 


■L( t ) I k - L(t)H(t) J 


(y^e 3R m , y 2 £ 1R *) which puts (2.1) in the form 


£(t) = (A u (t) - A 12 (t) L(t ) ) §(t) + g 1 (t) (2.2a) 

e f) ( t )= (A 22 (t) + eL( t ) A 12 (t)) rj( t ) + g 2 (t) (2.2b) 


with 

- g - 

CO 

= T~ X ( t ) 

>r 


T) 


_ Y 2 __ 


But this involves in solving two Riccati differential 
equations given by 


s L = A 22 L — A 2 ^ " £ L(A^^ ~ ^12^ (2.3a) 

• CO fsj 

~eH = HA 22 - A 12 + e HLA 12 ~e( A 11 -A 12 L)H . (2.3b) 



Again, since the slow state (or a linear combina- 
tion of and y£) gives no trouble (analytically 

or numerically), we need (for our purpose) only 
separate the slow state and allow the differential 
equation for the fast state y^ ^ 0 contain the slow 
state variables (see iviatheij [74]). To this end, we 
assume the following : 


(A~l). The matrices A^.(t), i,j = 1,2 are continu- 
ously differentiable. 

(A -2). There are constants ct. > 0 and c . . > O 

such that | |A ij (t) | | < c^ , 1 ( A ± j ( t ) | | < c?^ 

(A-3). There are constants c^ > O and c 4 > O 

such that any eigenvalue x (t) of A 22 ^^ 
satisfies : 

Re \(t) < - c 3 or Re x(t) 2: c 4 • 


Under these assumptions we will establish the 
existence and boundedness of a transformation 



which decouples the slow variable from the system 


( 2 . 1 ). 




Let E 


0 


» 



0 


then E' 


m 


0 


e I| 


Assuming that L is bounded T”^ admits of the fol- 
lowing representation. 

™I -eL~| 

„1 m 

T X = 

0 X k 



A 11 

■ A 12 




Denote A = 




f = 



< — 1 

Of 

< 

1 

A 22 


i 

! 

f 2 


The system (2.1) can be rewritten as 


E y = 

= Ay + f 

(2.5) 

<==> 

y = E“ 1 Ay + E” 1 f 


<==> 

ET" 1 ^ = ET" 1 E~ 1 ATT" 1 y + ET"' 1 E“ 1 f 


<==> 

_1 . _!• „l 

ET y - ET iT y 



= ET" 1 (E" 1 AT-T)T" 1 y + ET“ 1 E“ 1 f 


fl 

II 

V 

E^r( T "' 1 y) = ET” 1 (E“ 1 AT-T)(T" 1 y) + (ET” 

U v 

’ 1 E" 1 f) 


1 

(2.6) 

Substitute z = T y 

(2.7a) 


g = ET” 1 E” 1 f 

(2.7b) 

in (2 

.6) to obtain the equivalent system to 

(2.5) 


and hence to (2.1) : 



Az + g , A = ET“ 1 (E~ 1 AT - T) 


( 2 . 8 ) 


Now, 

El~ l E -1 


m 


0 


(2.9a) 


and 


A 


A 


•~ la 2 1 

eA ll L+A 12~ eL 

21 

e A 21 L -f A, 




22 


(2.9b) 


In the system (2.8), the slow variable is decoupled 
from the system, if and only if the upper-right 
block in A is zero? i.e., we put 

o 

sL = — L A 22 + A^ 2 + s(Aj^ ~ A 2 i)L (2. 10) 

Our task now is to show that for sufficiently small 
s , there exists a bounded continously differentiable 
matrix L(t) satisfying (2.10). We choose the 
boundary conditions for (2.10) in such a way that 
the boundary layers in its solution are eliminated. 
Thus we are justified in seeking a solution of (2.10) 
in (an asymptotic) series form. However, the eigen- 
values of A 22 (t), the system matrix of (2.10), are 
not grouped, i.e., with respect to their signs of 
real parts* which is necessary for the dichotomies 
involved. Thus before seeking an asymptotic solution 
to (2.10), we first order (regroup) the eigenvalues. 



* 

From (A-3) and the boundedness of A 22 , A 22 it 
follows that (see Coppel [20], Gingold [38]) there 
exists a non-singular continuously differentiable 
bounded matrix W(t) with bounded derivative such 


that 


W(t) A 22 (t) W~ X (t) = D ( t ) = 


where 


D x (t) 0 
0 D 2 (t) 


(2. 11a) 


ReX (D x (t)) < - c 3 < 0 

ReX(D 2 (t)) > c 4 > 0 . 


(2.11b) 

(2.11c) 


Regarding this, we let 


L = MW 


( 2 . 12 ) 


to reduce (2.10) to : 

e M = -MD + A 12 W“ 1 + e(A 11 -MWA 21 )M - e MVfv'T 1 (2.13) 
We seek a solution of (2.13) in the form : 

M(t) = 2 s^M . ( t ) + e R(t) (2.14) 

j=o 3 

In fact we are interested in an approximation 
0(e 2 ) , i.e. , N = 2. 

Substituting (2.14) in (2.13) and comparing 


like powers of e , we obtain : 



(2.15a) 


M 0 = A^WV 1 = A^A^W- 1 

M. = (~k + A,.M - M WA 01 M - M WW" 1 )D“ 1 (2.15b) 

1 v o 11 o o21o o 

and so on. 

In general, the remainder term R satisfies an 
equation of the form ; 

eR = - RD - M n _ 1 + f 1 (M 0> ...,M N-1 ) 

t £ f 2 (R » , • • • > __ » s ) * (2.16) 

Using (2.16), we can get an approximation of any 
order, i.e., 0(e N ) for any N. We are interested 
only in N = 2. In this case, 

f l = < A n “ M o WA 21 > M 1 " M l( WA 21 M o + ^VT 1 ) (2.17a) 

f 2 (R) = ( A n ~ j ^o WA 21 ” 6M 1 WA 21 )R 

" RWA 21^ M o +eM i + s2r ) ~ e RWW" 1 (2.17b) 

and (2.16) reduces to 


eR = - RD ~ M x + f 1 +ef 2 (2.18) 

The above discussion shows that our aim is 
accomplished provided that the solution R of (2.18) 
remains bounded. 


Let <P(t,s) be the matrizant of the system 


x = D^(t)x and ¥(t»s) 
system x = D2(t)x. 


be the matrizant of the 


CENTS U. U5RARV 

• \ '!• k T,. KANPUR 


4 r% a > * 



M 1 (t) + f x 


(2.19a) 


Let N (t) = - 

0 


N(t) = N 0 (t) + ef 2 . 


(2.19b) 


Define an integral operator S : IR — -> IR by 


t 


S<t) = ± / N(s) 


<P ( s , t ) 0 


0 


0 


ds 


-if n ( s ) 


0 


0 


o f(s,t) 


ds 


( 2 . 20 ) 


Lemma 2.1 


Let R = { R : j Jr| j < P] for some P > O. 
If (A-l) ~ (A»3) are satisfied then, there are 
constants K^, K^, K^, e ^ such that 


(a) | | N 0 (t) | | < K x , 

(b) | | N(t) || < K x + ePK 2 + ^^3 , 

(c) || f 2 (R x ) - f 2 (R 2 ) || < K 4 || ft r R 2 || for R 1 ,R^R 

and 0 < e < e ^ » 


(d) R £ R is a solution of (2.18) iff it 

satisfies the integral equation R = S. 


Proof 

(a) - (b) follow from the boundedness of 
Aij ( t ) 3nd ( t ) , M Q (t), M 1 (t) on [0,l] . (c) 



follows from using the matrix identity 


R 1 WA 21 R 1 ~ R 2 WA 21 R 2 

= (R^ - R 2^ A 21 R 1 + R 2^ A 21^ R 1" R 2^ 

since, ^2^ R 1^ " ^2^ R 2^ 

= (a x1 - m 0 wa 21 - 6 m 1 wa 21 )r 1 

- R 1 WA 21 (M q +eM x + e 2 R x ) ~ £ R i^ w 

“ ^ A ll ~ M o WA 21 ~ e M 1 WA 21 ) R 2 

- R 2 WA 21 (M q + eM x + e 2 R 2 ) ~ eR 2^ W ^ 

= (A u - M q WA 21 - eM 1 WA 21 )(R 1 ~ R 2 ) 

- (r 1 - r 2 )[wa 21 (m 0 + sm 1 ) + ewvr 1 ] 

- e 2 [R 1 WA 21 R 1 - R 2 WA 21 R 2 ] ’ 

and W as well as W* 1 are bounded. The ’if’ 
part in (d) follows from differentiation and the 
'only if’ part that from variation of constants 

formula. 

Lemma 2.2 

If (A-l) - (A~3) are satisfied then there 
exists e* > 0 such that the integral equation 
R = S has a unique solution in R, for 0 < e <e 

Progf 


By contraction mapping principle, the lemma 



will be proved if S happens to be a contraction 
mapping. 


Now, S(R 1 ) ~ S(R 2 ) 


1 1 r . <p(s,t), 0 j 

~ / [N + efp Ri) ds 

e o 0 21 q o 


-■/ [N 0 +ef 2 (R 1 )] 0 s(s>t) ds 

u 


“ / [N 0 + 6 f 2 (R 2 )] 


<P (s ,t ) 0 


+ 7 f IX + e f 2 (R 2 )3 


e t ° 


S (s,t) 


e / e C f 2^ R l^ " f 2^ R 2^ 

o 


<P (s ,t) 0 




f 2 (R 2 ^ 


* (s,tU 


since N Q does not depend upon R. Hence from 
lemma 2.1 and theorem 1.1 it follows that there 
exists e 2 > 0 such that, for 0 < e< e 2 ;<min(e o ,e^) 



| (StR.,.) - S(R 2 )| | 


< [ } K e ” c 3( t " s )/e ds + / K e" c 4^ t ~ s ^ e ds] 

o 1 0 


x K, 


R 1 “ R 2 


> e max(K,K ) .K. 

- in^rfp" • n*i -Ml • 


Hence 


| ! S(R X ) - S(R 2 ) || < ± | | Rj_ 

min (c 3 , c 4 ) 

2 K 4 max(K,K 0 ) 


for e < e. 


Ro II 


( 2 . 21 ) 


Also, from Lemma 2.1, || S 


< ~ [ / K e ” c 3^ t “ s ^ e ds + / K e ” c 4^ t '~ s ^ e ds] 


x (Rj_ EPKg + e 


max(K,K Q ) 

min(c 3 ,c 4 ) 


(K, + ePK 0 + e 2 P 2 Kj. 


(2.22) 


We choose 

2max(K,K 0 ).K 1 

p _ — -- — — 

min(c 3 , c 4 ) 



* min( Co , c , ) K 0 min(c„, c,) 

s = min( e , £, , ~ ♦ ' > • 

° 3K 2 max(K,K 0 ) 2K 4 .P 2K 4 max(K,K 0 ) 

Then, from (2.21 ) and (2.22) it follows that 
| jsj| < p and S is a contraction mapping on R. 

From the discussions above and the lemmas 2.1 
and 2.2, we obtain the following theorem : 

Theorem 2.1 

Assume that the system (2.1) satisfies the 
assumptions (A-l) - (A~3). Then there exist an 
e * > 0 and continuously differentiable matrices 
W( t ) ,M Q (t) ,M^(t) bounded on [0,l] such that 

W(t) A 22 (t)w~ 1 (t) = diag(D 1 (t), D 2 (t)) 

Re A (D x ) <~c 3 , Re \(D 2 ) > c 4 

and for 0 < s < e*, 

M(t) = M 0 (t) +eM 1 (t) + 0(e 2 ) 

the system (2.1) is transformed to 

Zf = [A lx (t) - M(t)W(t)A 21 (t)]z 1 + g x (t) 

e z 2 = A 2i (t)z 1 + [eA 21 (t)M(t)W(t) + A 22 (t)]z 2 + g 2 (t) 

= Zj_ + M(t)W(t)z 2 

^2 = z 2 



Lemma 3.1 


Consider the system 

ex(t) = [B ( t ) + e Y f(t,s)Jx(.t) (3.1) 

where j| T ( t ) j J c,_ s for every t £ [0,l] and 
y > 0. Let cp^tjs) be the matrizant of (3.1). 
Corresponding to this consider also the system 

ex(t) = S(t ) x (t) 0 < t < 1 (3.2) 

Let <P(t,s) be the matrizant of (3.2). Assume 
that | | B ( t ) | | , J | B ( t ) | | are bounded and 
Re X (B(t)) < -a < 0. Then there exist positive 
constants e*, oc^ and such that for all e , 

0 < e < e* , 

|| ^(tjS) - <F(t»s)|i < e Y K 5 e" G l (t " s)/e ' (3#3) 

for t >_ s. 

Due to the exponential dichotomy (Coppell [2l]) 
if the stability condition is replaced by 

Re A ( B ( t ) ) > a 2 > 0 , 

then, the inequality corresponding to (3.3) will be 
ll^t.s) - «'(t,s)| | < ^ e -M=-t)/e for t < s. 

V o « J 

Due to Lemma 3.1, it suffices to consider the fast 



system of the form 


ex = B(t)x + f(t) 0 < t < 1 (3.5) 

Assume that B(t) is conditionally stable and 

l|B(t)||, ||B- 1 (t)||, i|f(t)|| are bounded. Let 

{h } be a grid defined on [0,l], h = maxfh }, 
n n-1 n n 

t n = . Z h j > t 0 = °» t N = lf n = 

j=o J 

Assume that 

s 

min{h n } > max { y-gy , 2 e j E | } . (3.6) 


In fact, for standard methods to be accurate 
for the original system (2.1), one has to take, in 
general h << e , which is impracticable. Again this 
type of grid-generation with h<< e also affects the 
solution causing more round-off errors. Hence we 
are justified in assuming the minimum of the grid 
lengths to be greater than some number of 0(e). 
However (3.6) permits us to take very small grid 
lengths, since e<< 1. Let the numerical solution on 
the grids be denoted by u^ = u (t^) and similarly 
B(t n ) = B . In fact, in the course of error estima- 
tion we only need an assumption weaker than (3.6) , 


namely. 




< 



e 


(3.6* ) 



We will apply Euler's Implicit Scheme to solve sys- 
tems of the form (3.5). This particular method is 
chosen due to its unconditional contractivity proerty. 
This property holds when the one-sided Lipschitz con- 
stant v(t) < 0 (see Dekker et al. [22] s - v see also 
Kreiss et al. [68]). 

The difference equation corresponding to (3.5) 

is : 


h n ^ u n+l ~ u n^ 


B n+1 u n+l + f n+l 


or , 


e 


= (B 


n+1 


h n ^ u n+l + ^n+l 


(3.7) 


Due to the conditional stability, B n+ ^ is non- 
singular* thus (3.7) can be rewritten by 




B n+1^ 1 " h n B n+1^ u n+l + f n+l # 


In view of (3.6) or (3.6’), the matrix 


e i 

I - r— B is invertible. Hence the difference 
h„ n+l 
n 


scheme (3.7) is equivalent to : 


u 


n+l 


e , 

(X - 

B" 1 )~ 1 B 
h„ n+l' 

n 

n 

e 

- t 1 * h” 

O' 1 B n)l 


u 

n+l n 


n 


(3.8) 



But IU - — 

n 


= |l + 


+ & b :-> 2 + ••• 


n 


h n n+1' 


(see 


- 1 + I h n B n+1 I + I B n+1 ^ + *** 


(1 - I ff~ I)" 1 (due to (3.6')) 

n 


Hence (3.8) implies that 


u n+l « - h 


B 


-1 i 

n+1 ' 


" 1 " ^ 1 B nil I 


v + 




n+1 1 


1 ~ h n l B n+l 


1 i l f n+l^ 


-1 .-I 


Letting a n+1 = |3~ +1 


0 ( 1 ) 


0 = e/h n < 1 (<< 1) 


the above inequality reduces to 


u n+ll ^ 


0 


f 


a n+l ~ 01 n 


u + 


n+1 1 


a n+l “ 0 


a 


Let y 


n 


n+1 

0 


- 2 > 0 


(due to 6' or 6) 


0 < p 


f 


n+1 


n 


a 


n+1 


[54]) 


(3.9a) 

(3.9b) 

(3.9c) 

(3.10a) 


2© 


0 ( 1 ) 


(3.10b) 



Then 


1+Y 


n 


LA 


a , •* 

! + _n±i „ 2 

9 


JSL 


a 


n+1 


9 


a . , f , , 

/ n+1 r> \ n~rl 1 

~ ^ a n+! ~ 


a , , 

i + ~a±i - 2 


n+1 


I / 9 


If 


( a n+l ~ « )/e 


n+1 1 


a 


n+1 


0 


Hence (3,9c) can be written as : 


u n+l * ~ 1+Y 


n 


u i + LsA 

n 1 1+Y n 


(3.11) 


Applying theorem 6.2 in Chapter 1 we obtain 


n-1 , 

u„| < , 0. + ( % ttT~) ju. 


n 1 


max 

0 < 1 < n-1 K j 




(3.12) 


But |3j 






2B/ hj 


1 + y j 


e/n. 


B 


j+1 


e/h j 


Hence , 


u_ 


max 




1+1 


n 


0<j<n-l 


B 


-1 ,-l 
1 + 1 ' 


2e/ hj 



which shows that ju j is bounded. Similarly, the 
error-estimate can also be derived. Let the solu- 
tion x(t) be twice continuously differentiable, 
or, in other words, B(t) and f(t) are so. Let 
us denote x n = x(t n ), the exact solution at the 
grid points t ' s. By Taylor's theorem, 


h n ^ x n+l x n^ ^n+1 x n+1 + ^n+1 + T n+1 

with l T n+1 i < \ h n max |x(rj)J. 

t n <^<t n+ i 

Then e n = x n “ u n satisfies : 


e 

h^ ^ e n+l e n^ = ^n+1 e n+l + T n+1 


(3.14) 


Now, (3.14) is of the same form as (3.7); analogously 
it follows that 



< 


+ 


max 

0<j<n~l 



n~l e / h . 

3=o jBjij- 1 ~e/h. 


2e/hj 



max 


h i 


T j+1 


0<j<n-l 


R 


-1 |-1 



Writing e = max j e n | , we obtain. 


e 


0(h 2 ) , if h | 1 3 | ) >> 1 
0(h) , otherwise 


We summarise the above discussion in the follow- 
ing theorem. 


Theorem 3,1 

Let B(t), f(t) be twice continuously differen- 
tiable j j 1 B ( t ) | j , ||f(t)|| are bounded. Let h n 
be a grid satisfying (3.6) (or (3.6')). Then the 
difference scheme (3.7) solves (3.5) with u bounded 
and error 

e = 0(h 2 ) if h | j B | | >> 1 ; otherwise e = 0(h). 


Remark 3.1 

If instead of Euler's Implicit scheme, we apply 
Euler's Explicit scheme and integrate backward, i.e., 
use the scheme 


(B 


n 




h n u n+l ^n 


(3.16) 


then as previously, we get the estimate (analogous 
to (3.13)) 



u 


i\ T ~n 1 


< max 


L N-j 


°- <Jln - 1 iB^jl' 1 +2^ 


n-1 

+ TC 

3=0 


s/h 


N-j 


b ~ 1 . 

i'!~3 


+ s/h 


u N j (3.17) 


N-j 


and the corresponding error estimate (analogous to 
(3.15)) : 


el < max 
e N-n 1 0< j<n-l 


h N~j T N~j~l 

1 + 2e 


n~l 

+ TC 


'N 1 


^ ° 1 + e 


(3.18) 


Again, here e = 0(h 2 ) if h|}Bj| >> 1, 
otherwise e = 0(h). 


Remark 3.2 

A look at Theorem 1.1 reveals the fact that if 
Re X (B(t)) < 0, the forward integration is stable, 
whereas, if Re X(B(t)) > G, then the stability 
holds in the negative direction? suggesting backward 
integration. Due to this reason, any numerical 
method which tries the solution of (3.1) should 



integrate forward or backward depending upon the sign 
of the real parts of the eigenvalues of the system 
matrix B(t). 

The above error estimates and Remark 3.2 suggest 
that we employ the Euler's implicit method and inte- 
grate forward when the eigenvalues of B(t) have 
negative real parts and employ Euler's Explicit method 
backward when the eigenvalues of B(t) have positive 
real parts. If B(t) have eigenvalues some with 
positive and some with negative real parts, a trans- 
formation (as in equation (2.11)) of 3(t) must be 
done to separate the eigenvalues with negative or 
positive real parts before applying the methods. 

Note that for such a transformation it is necessary 
to assume that on the whole of the interval [0,1], 
the spectrum can be divided, i.e., any eigenvalue 
X (t) of B(t) must be either having real part 
negative or positive throughout the interval. In 
practice, very often the contrary happens, i.e#, 

Re ( X (t)) changes sign, when t varies through 
[0,1]. To overcome this, we assume a weaker state- 
ment to hold. Let { b i } be a partition of [0,l] : 

0 = b„ < b, < ... < b, = 1. 

O 1 K 

Denote — |Jb^, b^_^), i — 1 , • « • , k—2 , ,1 = • 



(A-4) Assume that on each sub-interval the . 

eigenvalue criterion, i.e*, any eigenvalue 
X(t) of B(t) satisfies Re \(t) <_ -c^ or 

Re \(t) > c 4 

for some positive constants c^ and c^ 
is satisfied. 

By (a- 4) we are permitted to represent 
S(t) = diag ( B 1 ( t ) , B 2 (t)) (3.19) 

where the block B.(t) have eigenvalues 
with real parts of the sign (-1) through- 
out the sub-interval {3^ , i = 0,1,..., k-1. 
That means on each of the sub-intervals 
fA we must get transformations w i (t) 
which transforms B(t) to the form as 
in (3.19). 



2.4 Decoupling the Difference Equations 


In the preceding sections we have seen that 
Euler's Implicit method is 0(h) or 0(h) accurate 
for the fast variable y 2 « The slow variable z ^ can 
also be solved by Euler's Implicit method provided 
that we take h n smaller than the Lipschitz constant. 
This suggests to apply the scheme first and then de-- 
couple the variables if necessary. 

In this section we study this approach of 
taking the difference equations and the possible de- 
coupling transformat ion( s ) , applied on the difference 
equations . 

So, consider the system 


Ex = A(t)x + f(t) 0 < t < 1 



1 

M 

O 

1 


f 1 

with E = 

m 

f = 



„ 0 eI k„ 

i 

? 

i 

f 2 


x 


X 


1 



x 1 em m , 


0JR 



(4.1) 


A(t) 


A 11 A 1 ^ 


,21 .22 

A A 


similarly 



Assume that (A~l) of § 2.2 holds. The discretiza- 
tion by Euler's Implicit scheme yields : 



1 

x n+l 

- x 1 
n 


Fa 11 

n+l 

A 1 !," 

n+l 


1 

x n+l 


,1 

” f n+l 

E 

x 2 

n+l 

- x 2 
n 


A 21 

n+l 

22 

A”, 

n+l 


x 2 

j x n+l 

+ 

f 2 

n+l 



— i 


— 



L_ - 




(4.2) 


(4.2) can also be written in the abbreviated form as 


h; E " A n+l> x n+l 


r* Ex + f , 
h n n n+l 


(4.3) 


Introduce a sequence of transformations ( T n } by 


defining (formally) 


n 


m 


eM. 


n 


0 


(4.4a) 


whose formal inverse is 


As previously, put 



'm 


0 


-eM 


n 


(4.4b) 


(4.5) 


to obtain (from (4.3)) 



E z n+l 

h T n+l^ h n A n+1 ' 

- S > _1 £ T n 2 n 


” h n £ T n+l^ h n 

A n+1 “ f n+l 

But h n 

A .. . c — 

n+l 

h A^ 
n A n+1 

h A 21 

L n A n+1 

h A 12 

-I n n n+l 

m 

h a 22 

h n rt n+l ~ 1 


(4.6) 


el, 


which implies that 


< h n A n+1 - 


-Do 


-C, 


(4.7) 


where , 

c 3 = (h r 
Dyi = (h n A 


c 4 ~ A 


D, 


11 

'n+l 

- z m - h n An+N+An+t * 

sl k>“\ A nli rl 

22 

n+l 

- eI k - h A n+l< h n Cl ■ 

T >>~1 K A 12 N-l 
^m h n A n+1 

11 

'n+l 

" I rd \ A n+l °4 


r — 1 

CM + 
CM C 

- 6l k )_1 h n Cl C 3 • 



Note that if the inverses above exist, then they 

are bounded. The inverses exist provided that 

min h |A^* | < 1. From (4.6) and (4.7), we find that 
n 



1 



— 



z n+l 


~°3 ~ e M n+1 D 3 

- e \+i- 


i 

z n 




* D 3 M n + eC 4 + 






+ ^“n+l U 4 



„ z 2 

Li n+1 ~.i 


i_ s °3 

e2 °3 M n - 

•sw m 


2 

z n 


+ 


(4.8) 


C 3 + sM n+l D 3 


e D. 


0 A - sM ..D. 
4 n+1 4 


D 


4 


Ti+1 


n+1 


To claim that (4.8) is a (partial) decoupled form of 
(4.6), we must put the upper-right block zero, i.e., 

eM n+l (D 4 " D 3 M n ) ” c 3 M n + C 4 = 0. (4.9) 

This is a well-known Riccati algebraic equation. If 
we take M ^ as the unknown, consequently, the 
Euler Implicit is applied forward, then (4.9) is an 
ill-conditioned equation, clue to the appearance of e 
with M Instead, if we integrate the system with 

Euler's backward implicit scheme, then equation (4.9) 
might be interpreted as an equation for the unknown 
matrix M n and (4.9) is no more ill-conditioned • 
moreover M n remains bounded. 



Matheij [74] considered a class of multistep 
methods to the system (4.1) with the stability 
criterion that Re A(A 22 (t)) < 0, for all t £ [0,l],‘ 
he introduced such a sequence of transformations (4.4). 
The above interpretation of (4.9) reveals that one 
must use Euler's backward difference scheme, which 
is a member of the class of methods used in [74 j. 

In the following we study the corresponding 
decoupling of the difference equations for the more 
usual decoupling (Lyapunov) transformations used, 
for example in Kokotovic et al. [62] (cf. § 2 . 2 ) 

This approach was also taken by O'iLalley and 
Anderson [86] but for the constant coefficient case, 
i.e,, when A(t) does not depend upon it* 

For this purpose, we take the difference approx- 
imation (4.3) and define a sequence { S n } of trans- 
formations by 


where L n and H n are two matrices of dimensions 
kxm and mxk respectively to be defined (deter- 
mined) later. 


m 


e H 


n 


-L I. - eL H 

n k n n 


(4.10a) 



Th.e formal inverse of $ n is : 


.-1 

) 

n 


T m - H n L n 


n 


s H 


n 


(4.10b) 


Using the transformation (4.6) we put the difference 
equations (4.3) in the form below. 


Let z n = S n x n 


Then (4.3) admits the representation 
Ez n-U " ( ES W lE V z n + < ES n+l E > f n + l 


where B = (E - h n A n+1 )' 


~ C 3 C 4 


D, 


-D, 


'3 4 

with » C^, as used in (4.7). 

Using (4.13) in (4.12), we obtain 


1 

z n+l 


r c 

l i 

F 2 ~ 


'^1 

2 

eZ n+l 


F i 

C 2 


! 

N 

D fO 
l 


^ ” I m +H n+l L n+l ^ C 3“ e H n+1 D 3 


(4.11) 


(4.12) 

(4.13) 


sL n+l C 3 +eD 3 


^ I m" H n+l L n+l' )C 4 +eH n+l D 4 


e L n+l D 4" eJ 4 


n+1 


'n+1 


(4.14) 



where , 


” S H n+l L n+l^ C 3 e H n+1 D 3 
e (V EH n + l L n + h C 4 L n +E2H n + l D 4 L n 


c„ = 


e2( - L n + l C 3 +D 3) H n +e2 ( L n + l C 4- D 4>( I k- eL n H n ) 


eL n+l C 3 +eU 3' 


’^ L n+l C 4~ D 4^ L n 


8 ( I ni" eh n+l L n+l) C 3 H n“ e2H n+l D 3 H n 
+ £ ( I m~ eH n+l L n+l) C 4( I k'' 6L n H n) 


+ e !i x ,D,(L- eL H ). 
n+1 4 V k n n' 


Since the transformations of the type (4.6) contain 
two unknown matrices L and H, we can afford to de- 
couple (4.3) completely. It is also necessary to 
do that since we need determining L as well as H. 

In (4.14) it is obvious that, for this purpose, we 
must put = 0 = F Again F 1 = 0 simplifies the 

equation F 2 = 0 as well as the diagonal blocks 
and C 2 . Hence, if L n , L n+1 and H n , H n+1 satisfy 
the equations : 

L n+l (C 3 +€C 4 L n ) ~ ^ D 4 L n +D 3 ) = 0 (4.15a) 

H n+l^ D 3~ L n+l C 3^ “ C 3 L n “ C 4 L n') H n L n = 0 ( 4 * 15b ) 

then the decoupled equations are : 



(4.16a) 


z n+l 

e2 n+l 

where 


-( c 3+ eC 4 L n )z 


n 


(~D 4 + e L n C 4 )z 


n 


Q 


n+1 

.2 

J n+1 


is 


+ ^n+1 

+ g^ +1 (4.16b) 

the non-homogeneous term in (4.14). 


Equation (4.15a) is a well known Riccati 
Algebraic equation which can be solved for L n+1 
and leaves bounded. Note that we start the 

recursion with L q = 0. (iviatheij [74]). The same 
remarks apply to equation (4.15b), provided that we 
have already solved (4.15a). The recursion (4.15b) 
might also begin with H q = 0. 

But if we take the Euler’s implicit scheme 
backward, i.e., with the same scheme we want to 
integrate in backward direction, then it follows 
that we must solve the equations (4.15) for and 

H . In that case, the matrix equations (4.15) become 
ill-conditioned. The solutions, thus can not be 
expected, to be bounded. 


2.5 Comments on Decouplings and Difference Schemes 


In the preceding sections, we have seen that 
Lyapunov type transformation decouples the slow and 
fast variables completely, whereas Matheij type 
transformation decouples the slow variable from the 
system. For numerical purposes, however, the latter 
has advantages, since it needs only one Riccati 
Differential Equation to be solved, whereas the 
former needs two. This also holds for the approach 
taken in Section 2.4, where a difference scheme is 
first applied and the decoupling next. In this 
approach, Euler's Implicit scheme works well for 
the former transformation (a sequence of trans- 

formations) and Euler’s Explicit scheme is appropri- 
ate for the latter (Matheij [74]). It is also imp- 
ortant to know that in which direction the integra- 
tion is to be carried out, backward or forward. Due 
to these restrictions, the methods have limited 
scopes, for example, they are appropriate for initial 
value problems. The direction of stability can also 
not be neglected (Theorem l.l). 

Therefore, we look at decoupling the states 
before application of a difference scheme. Decoupling 
might be done by any of the transformations ; 



Lyapunov type or the one introduced in g 2.2. However 
we use the later for obvious computational convenience 
Here again we need some clarifications as to how to 
solve the decoupled systems for given boundary condi- 
tions. If the boundary conditions are separated, 
i.e., they are given in separate sets for the fast 
and slow variables, then one might solve the decoupled 
slow variable first (e.g. by a centred difference 
scheme), and then use the results for the computation 
of the fast variable. If this is not the case, then 
one faces trouble as to generating an appropriate 
grid. For the slow variable, a grid $ h n ) is 
required, where h is less as compared to the 
Lipschitz constant of the slow system. On the other 
hand one might like to have h to be greater compared 
to the reciprocal of the norm of the system matrix of 
the fast system (see Theorem 3.1) to guarantee 
e = 0(h). This situation can be remedied by using 
a stretching (not necessarily uniform on the whole of 
[0,l]) of the independent variable. In fact, we need 
the stretching where the spectrum-structure changes, 
(see the Hybrid algorithm in $3.2). 

We have also remarked that the assumption that 
the eigenvalues with positive or negative real parts 
of the fast system need not be assumed to be separated 



on the whole of [0,l]. (see (A~2)). Computationally, 
then we must determine the subintervals [b^, b^ + -^] 
on which this happens, i.e., the structure of the 
spectrum does not change within any of the intervals 
(b i> hi+l 5 - To this end, we select a partition 

b i » 0 = b o < b l •••<•••< b b = 1; find the 

1 1 

transformation W(bt) and block-diagonalize the 
system matrix of the fast system. This can be 
achieved by the usual QR-algorithm, or/and its modi- 
fication by Kreiss et al. [68], We check, whether 
the spectrum structure does change for different 
b^'s. If it changes for b^ and b b + ^ > we might 
need another point in [b^, b i+iJ* that, if 

the spectrum structure is not uniform throughout [0,l], 
then there are points b^’s, where the change of 
the structure is smooth, and our aim is to determine 
such points, (see §3.2 for details). 

Our last remark concerns the global error. For 
the separated sets of boundary conditions, the global 
error can be estimated (see Dekker et al. [22]). But 
for the general boundary conditions, the global error, 
as such can not be estimated beforehand. Because, 
the structure of the spectrum is not known beforehand, 
and thus , the difference method is not known completely 
beforehand (except conditionally), condition number 



of the error matrix can not be determined, and con- 
sequently the global error can't be estimated (see 
Kreiss et al. [68]). 


Chapter 3 


NUMERICAL DEMONSTRATIONS 

3.1 Introduction 

In the previous chapter we have analyzed, the 
decoupling transformations and the difference 
scheme(s) for linear singularly perturbed systems. 

We have described the grid-generation and the 
numerical method that combines decoupling and dif- 
ference schemes. The error analysis shows that the 

p 

method is accurate upto 0(h) (or 0(h)) at least 
locally. The global error analysis is not included, 
since, in general, it is not possible at this stage 
to do anything conclusively. 

We realize that though a 'good theory' is a 
must, yet it does not conclude once and for all that 
a numerical method would never be unsuccessful. One 
loophole is the roundoff errors, which we do not 
include in the error analysis. Round-off errors 
depend heavily on the complexity of a method and the 
number of functional evaluations necessary for work- 
ing it out. That is why it is instructive to do 
some numerical experimentations with any proposed 


method . 



This chapter aims at fulfilling this gap bet- 
ween theory and practice. In § 3.2, the hybrid 

algorithm is presented, where we summarize the 
hybrid method in an algorithmic form to make imple- 
mentation easier. In s§ 3.3, we take some examples 
to demonstrate the effectiveness of the method. The 
last section comments on the usefulness, effective- 
ness and the limitations of the algorithm. 


3.2 Hybrid Algorithm 

To facilitate understanding of the algorithm 
let us write down the problem we are considering 
and the underlying assumptions once again. We con- 
sider the linear system of singularly perturbed 
two-point boundary value problem : 

+ f x (2.1a) 

ey 2 = A 21 y l + A 22 Y 2 + f 2 (2.1b) 

f or 0 < t < 1 with the boundary conditions 


n = A nn + A i2>2 


Yl< 0) 

+ B , 

y x (D 


1 

y ; ,( 1) 


(2.1c) 


where 


i t 


denotes d/dt, t - often called time, 



e is a small (compared to l) positive parameter, 


A^j(t,e) are smooth real-valued matrix functions 

of dimensions n. x n. , fls are smooth real valued 

^ J J 

vector functions of dimensions n.. where n. = m, 

3 J 

if j = 1 and n j = k, if j =2; B 0> are time- 

independent real matrices of dimensions (m+k) x (m+k); 

is an m+k dimensional time-independent real 
vector so that (2.1c) represents m+k linearly inde- 
pendent boundary conditions. We assume that : 

(AL~l) There exists a partition . {c^} , 

0 = c Q < c-^ < ... < Cj = 1 such that on 

each of the sub-intervals [ c^, c^ + ^ ] there 
exists a constant > 0, so that the spec- 

trum of A 22 ( t , s ) , i = 0,..,J~1, can be 
represented by the union S_U S + where 

Re A < _ K i , if ^ £ S_ and Re A > K ± , if 

X e s + , 

(AL-2) There are constants > 0 such 

that | i A ± j | | < K i j and Ma^M < , 

1 9 3 = 1 » 2 . 

Note that the smoothness of A^j and fj 
are included in the formulation of the 
problem, though, of course, we only need them 
to be twice continuously differentiable. 



An outline of the algorithm follows. 


Our first task is to realize (AL-l) numerically 
i.e., to obtain the partition { } and transform 

A ^2 to the block form D (see (2.11) of Ch. 2) on 
each of the sub-intervals [c^, c^ + ^]. We do it 
recursively. Given a partition point c^, the next 
partition point c j_ + j_ is determined as follows : 

(A) If c • , j = 0,1,..., i have been previously 
determined to be partition points, then com- 
pute the eigenvalues ^j(°^+), where 
Aj(t) =Xj(A 22 (t)), using say, QR-algorithm, 
thereby forming new sets of eigenvalues 
S_ (c^+), S + (c^+) in addition to the former 

ones S (c- ) and S (c.). 

- i +1 


Note that S(c-;+) means lim S(c.+h). 

h*0+ 1 

Numerically c^+ is to be interpreted as 
some point greater than c^ but sufficiently 
close to it. (Typically we take c^ = c^+e/10). 


(B) 


Let c = c^ + a be 

c i+l J w ^ ere a = - 

take some trial value, 


a trial value for 
c.^ (if i = O, we 
typically = l/l6). 


Compute the eigenvalues (c) and determine 
the sets S_(c) and S + (c). 



(i) if no sets S change (compared to S (c^ ), 
and S + (c i )), go to ,(iv). 

(ii) if the sets S change, mark c as a possible 
partition point. Before going to the next 
step, it is worthwhile to check, whether 
there is a partition point in between c^ 

and c , 

(iii) try (i), ( ii ) for a = (c. -» C ^_^)/V 2. (It 
is advisable to try this again and again by 
taking smaller a's). It might also happen 
that we took a unnecessarily, i.e., the 
same structure of the spectrum is valid for 

a larger sub-interval. 

(iv) if (iii) does not add anything new, try (i), 

(ii) for oc = 2(c^ ~ c^_^). Go to ( C) . 

The next task is to decouple the system 
in each of the sub-intervals [c^, c^ + ^]. 

(C) Compute the matrices W(t) and D(t) using 
QR~algorithm (see Wilkinson [115] also 
Kreiss et al. i.e., lemmata : 9. 1,9.2, [68] ) , 
on each of the sub-intervals [ , c i+1 ] , i.e., 

on the grid-points, (see Remark 2.1). 



Remark 2.1 


We first take a uniform grid in each of the 

sub-intervals [c^, c j_ + j_] (typically grid length 
= l/l6). Then determine the matrices W(t) and 

D(t) on these grid-points; Note also that the 

grid should satisfy (3.6') of Ch. 2. Compute : 

M o 

“ A 12 A 22 W 

(2.2a) 

ivi. 

== (A 1X M — M q — ~ MqVvI/V )D 

(2.2b) 

M = 

i\/i + e M, 

o 1 

(2.2c) 

B 11 

= A X1 - MWA 21 

(2. 2d) 

B 22 

= A 22 + eA 21 MW 

(2.2e) 

i — 1 

= f x - MWf 2 

(2.2f) 

Then 

the 

the system (2.1a - 2.1b) is transformed into 

form : 

*1 

= B ll z l + 9l 

(2.3a) 

• 

e 

~ A 21 Z 1 + B 22 z 2 + f 2 

(2.3b) 

with 

Yl = Z 1 + MWz 2 > z 2 = y 2 

(2.3c) 

on each of the sub-intervals [c^, 

Our task 


now, is to solve the system (2.3a - 2.3b) with (2.1c). 



Mote that the boundary conditions (2.1c) should 
also be transformed as well. This is done using 
(2.3c). Write the transformed boundary conditions 
as ; 

*i(o; 2^(1) 

+ C 1 = C 2 . (2.4) 

z 2 (0) z 2 (l) 

Note that the system (( 2. 3a) -(2.3b)) is an m+Ic. 
dimensional system, with (2.3a) being of dimension 
m and (2.3b) can be represented as two systems 


with dimensions kj., 

k 2i ^ k li + 1 

k 2i 

= k) on each 

of the 

sub-intervals 

^ c i * C i+1^ » 

•H 
* — ! 
M 

being the 

order 

of D, and 

k 2i * ‘ tha 't 

of 

D 2* 


We intend to solve the slow system (2.3a) by 
Trapezoidal rule and the fast system (2.3b) by the 
Euler's (backward and forward) schemes ". 

Hence according to the above paragraph, trapezoidal 
rule is applied for the first m dimensional system, 
Euler's forward Implicit scheme for the next k-^ 
dimensional system and Euler's backward Explicit 
scheme (see § 2.5) for the last k^^ 

system. This is done in (D). 



dimensional 



Remark 2.2 


For the methods to be accurate, the grid must 
;d accurately, i.e., tl 
satisfy s (see (3.6') of Ch. 2) 


be generated accurately, i.e., the grid { h n } must 


2s. 

h n 


n 


B, 


< 


_n 

s 


< 


3 


11 ' 


for all n. 


If possible, the condition that ^ n f^22^n+l^ ^ ^ 
also be verified. A sub-interval where the grids 
are unusually refined compared to other parts will 
be referred to as a stretching sub-interval, and 
its end-points as stretching points. 


(D) Denote 


6 . . 

ij 


1 

2 ’ 


3 — 1 ,2, , . . ,m 
3 = m+1, . . . ,mtk li 
j = m+k-^^+1, . . . ,m+k 


(The first subscript signifies the choice 
of 6.^. is dependent upon the sub-intervals 

~ c i ,c i+l*^ * ^ 

Denote A i = diag (5 i j ) . The (hybrid) dif- 


ference scheme is : 



u , . - u 
n+1 n 


h 


- n : A i B n u n + ^ A i^ B n+l u n+l 

+ A.G n + (I - A i )G n+1 

for t n 6 [ c iS c i+1 ]- (2.5a) 

where, we have abbreviated (2.3a - 2.3b) by taking 


B = 


3 

w 

b- J 
h- 1 

0 


1 

* — I 
O' 

1 

-1, 

~i 0 

9 G = 


r — 1 

CM 

< 

CO 

f 

e b 22 


f 2 


(2.5b) 


and the usual notation x n means (x : u, B, G etc.) 

x(t n ) , t n = h n * n = 0, . . . , N 

1 = 0 

The boundary conditions are added with the difference 
system as usual by taking 


z j_( a) 


22 (b) 


» U N ~ 


z 2 (a) 


_ Z 2 (b) 


(for details, see Keller [55]). 

(E) The numerical solution at the points not 

included in the grid is determined by using 
Richardson's extrapolation [ l]. 

Remark 2.3 

Note that stretching of the boundary layers 



etc. is done automatically in the form of grid- 
generation due to the conditions in Remark 2.2 
that are to be monitored in the grid-generation 
process . 

Remark 2.4 

An alternative approach is to solve Equation 
2.10 of Ch. 2 approximately by using asymptotic 
series expansions. For example 0(e) approxima- 
tion would look like ; L = L q + e + O(e^) 
with L q = A 12 A 22 

L 1 = ^ A ll “ L o A 21^ L o " L 0 ^ A 22* This a PP roxi " 
mation of L lets the system matrix of the decoupled 

equations look like 


A 


11 


LA 


21 


A 


21 


0 


A 21 L + A 22 


Note that this is to be done on the whole of [a,b] 
at once (contrary to the previous approach, where 
W ( t ) was involved.] Then we find out the trans- 
formations W(t) on the sub-intervals 

This approach is justified due to Lemma 3.1 in 


Chapter 2. 



3.3 Examples 


We have seen the difficulties associated with 
linear systems of singularly perturbed two-point 
boundary -value problems in § 1.1. The difficulties 

were, in general, two fold from the vi ew point of 
computations. As was noted there, computation of 
the solution should take care of the situations 
with large eigenvalues and situations with the 
largely varying eigenvalues. Again, asymptotic 
analysis of these problems reveals the fact that, 
very often, in the practical problems, there occur 
boundary layers near the ends t = 0 and t = 1 of 
the interval of integration [0,l], 

The corner layers or the interior layers are 
also observed to occur in the middle of the 
interval . 

In this section, we consider some examples 
where we meet one or many of the aforesaid situations. 
Figures following the examples show the numerical 
solutions. They intend to show that in the choosen 
grid the solution do not contain any oscillatory 
elements as well as they give ideas of how a solu- 
tion looks like. Only for modest e 's, these 
figures are drawn. The tables following the examples 



provide the absolute value of the difference between 
the exact/asymptotic solution and the numerical 
solution obtained by the hybrid method. In the 
tables the following abbreviations are used. 


SP 

ss 


NG 


Mad 


Mid 


Stretching points (end points of SS) 

Stretching sub-intervals, where the grids 
are made finer. 

Number of grid-divisions used in the cor- 
responding stretching sub-interval. 

Maximum of the absolute values of the 
differences between the exact/asymptotic 
solution and the numerical solution. 

Minimum of the absolute values of the 
differences between the exact/asymptotic 
solution and the numerical solution. 


For each of the problems, we give two tables; where 
the first table gives the choice of stretching points 
and the second, showing the (comparison of) results 
in terms of Mad. and Mid. 


From amongst many examples experimented, the 
following model example is choosen from O'Malley [ 82 ] 
for illustration purposes. The parameters a, p 
in the following example, chosen appropriately, give 



rise to different nature of the solution. We will 
also come across the difficulties associated with 
the eigenvalues of the system matrix, when the 
different cases come into picture. 

The model linear problem, with the constants 
(to be fixed later ) a, p is : 

sy + 2aty - apy = 0 - 1 < t < 1 (3.1a) 


with the prescribed boundary conditions given at 
t = + 1. (3.1b) 

The general solution of (3.1) is 

2 A 

y(t) = exp( - Si-)[c, D , 1 B (— ) 2 t) 
e x -x - 2 ^ e 

1 

+ c 9 D, (i(-^) 2 t)] (3.2) 

¥ e 


for arbitrary constants c-^, The constants c^, 

c 9 might be determined from the boundary conditions, 
i . e . , from 


y(±l) = exp(- &-)[c 1 D 


2s 


'1 — TTp 


1 . 

2* 


((+ f -) 2 


+ c D 1 (+ i(' ? ~) 2 )j. 
oP e 


(3.3) 


The asymptotic approximation of D n 1 s are : 

7 2 n 

D n (z) = exp(-~)z n (l + o(l) 

as |zj - <*>, for n = 0,1,2,... 



which might be needed in solving (3.3) for c ? . 

The asymptotic behaviour of (3.2) differ according 
to the values of the constants a and j3. We have 
the following four cases accordingly. 

Case 1 

a > 0, p ^ ~2n , n = 1,2,... . 

In this case the solution y(t) in (3.2) can be 
represented (asymptotically) as : 
for -1 <_ t < O, 

1 £ 

y(t) = (~t) 2 [ y(-l) - (-l) 2 y(l) + 0 ( 1 )] 

£ 

+ ((-l) 2 y(l) + o(l) 

for 0 < t < 1 , 

n _i _£) £ 

y ( t ) ~ — (y(-i) - (~1) 2 y(l) + °(l)) 


1 



tl + p/2 

k £ 

+ ( (-l) 2 y(l) + o ( l ))( — t ) 2 

and y (0) = 0(e^ 4 ) ; F (.) being the gamma function. 
This represents non-uniform behaviour at t = 0»an 
interior layer near t = 0. 



Case 2 


c: < 0, p ^ 2m, m = 0,1,2, . . . 

The solution is s 
for "1 < t < 0 , 

y(t) = — (y(-l ) + 0 ( 1 )) exp( a( l~t 2 )/e ) + 

t l+p/2 

+ 0 (e“U-«)^ ) 

for any 5 > 0. 

for 0 < t < 1 , 

-u+|) , 

y(t) = t ' (y(l) + o(l)) exp(a(l~t z )/ . 

and y(0) = 0(e a ^~ 

for any > 0. 

This shows the non-uniform behaviour of the solution 
at both of the end-points t = + 1 . 

Case 3 

a <0, p = 2m, m = 0,1,2,,.. 

The solution in this case is : 
for te [-1,0) U (0,1], 

y(t) =|[(y(l) - (-l)”Vi-l)) 

+ (y(l) + ( -l) m y (— 1 ) ) t m ] + o(l) 



Case 2. 


a < Q, p £ 0,2,4 , . . . 

Limiting solution as e-» 0 vanishes 
on - 1 < t < 1. 

Case 3. a < 0, p = 0,2,4, ... 

Limiting solution as e~ 0 is t 

+ y(-l)(-t)' 6/2 J on -1 < t < 1. 

Case 4. a > 0, £ = -2, -4, ... 

Solution becomes exponentially large 
on (-1,1). 

For computational purposes, we choose a, p, 
y(~l) and y(l) appropriate to each of the above 
four cases. In the following, the Example l.m 
correspond to the case - m. 

Note that a problem is always transformed to a 
system in the standard way. 

Example 1.1 a = 1, p = 2. (Case l) 

Explicitly, the problem is i 
ey + 2ty - 2y = 0 -1 < t < 1 

y(-i) = -l, y( i) = 2. 

Example 1.2 .« = -1, p = 1. (Case 2) 

1 • 6 • y 

0 

ey - 2ty + y = 0 
y(-l) = -1, y(l) = 2. 


1 < t < 1 



Example 1.3 


a - - 1 , 


p = 0 


(Case 3) 


i .e. , 

e y* - 2ty = 0 , -1 < t < 1 

y( -1) = -1, y ( 1 ) = 2. 

Example 1.4 a = 1, p = -2 (Case 4) 

i.e. , 

e y + 2ty + 2y = 0 -1 < t < 1 

y(--l) = --1, y ( 1 ) =2 

Note that the solution to Example 1.4 has exponential 
growth and damping in the left and right boundary 
layers respectively* and it is unbounded elsewhere. 

Due to this reason, we do not have graphical solu- 
tion for Example 1.4. Figures (3.1 - 3.3) show the 
numerical solutions for the Examples (1.1) ~ (1.3) 
with e = 10" 4 . Table ( 3.1 - 3.4) summarize the 
results of comparison of the numerical solutions 
with the asymptotic solutions. 

Note that the interval of integration is [-1,1] 
and not [0,l] which we adopted in discussing the 
hybrid method. This is handled by first transform- 
ing the interval [-1,1] to [0,l] before the application 
of the hybrid method. Alternatively, in the algo- 
rithm itself, any interval of the type [a,b] can 



be retained asking the user to provide the values 
[a,b|, and translating this[0,l] in the whole of 
the hybrid method into [a ,bj . 

To demonstrate the effectiveness of the hybrid 
method for higher -dimensional problems, the follow- 
ing example is chosen. This arose from the eigen- 
value problem in vibration of strings (Lord 
Rayleigh [71]). However, we will fix the eigenvalue 

X = tu 2 + 4btc 2 (+ 0(e 2 ) ) • 

in accordance with O'Malley [82, § 3.3 ]. 

Example 2 

Consider the two-point boundary-value problem : 

4 2 

~ d y d y 

2 . - — r = X y 0 < t < 1 (3.4a) 

dt* dt 2 

with the boundary conditions 

y(0) = y(O) = y ( l) = y(l) = 0 (3.4b) 


where , 

X = n 2 + 4 e-rc 2 . . (3.4c) 

The asymptotic solution for (3.4) is given by : 



y(t,e) = V2 n [ + e[ Sin-jrt 

u 71 71 

- (l-2t) cos Tit + e~ t//s ~ e -U-t)/e } + 0(e 2 )]. 

(3.5) 

The numerical solution for (3.4) with e= lCf & is 
shown in Figure 3.4 and the results of comparison 
of the asymptotic solution (3.5) with the numerical 
solution obtained by the hybrid method are summarized 
in Table 3.5. 

The experiences with the numerical experimenta- 
tion are presented in § 3 .4. 



Table 3.1a 


Choice of Stretching Points for Example 1.1 


£ 

SP 

b o 

b l 

b 2 

b 3 

1E- 

-2 

-1 

-0.05 

0.05 

1 

1E- 

-4 

-1 

-0.003 

0.004 

1 

1E~ 

■6 

-1 

-0.0018 

0.0008 

1 





T able 

3.1b 


Summary 

of the 

Comparison 

Results for 

Example 1.1 

0 


SS 

[b 0 > b J 

rb 1 ,b 2 ] 

[b 2 ,b 3 ] 



IMG 

12 

20 

8 

04 

1 

• — i 


Mad 

24E-3 

88E-3 

28E-4 



Mid 

OE-1 

IE-6 

OE-l 



IMG 

16 

16 

12 

IE -4 


Mad 

22E-5 

7E-5 

4E-4 



Mid 

IE -7 

IE -8 

OE-l 


NG 

16 

20 

12 

Mad 

2E-5 

22E-6 

IE -5 

Mid 

OE-l 

OE-l 

t 

rj 

o 


IE-6 




Table 3.2a 


Choice of Stretching Points for Example 1.2 


fj 

SP 

b o 

b l 

b 2 

b 3 

IE- 

-2 

~1 

-0.96 

0.964 

1 

IE- 

-4 

-1 

-0.992 

0.993 

1 

IE- 

-6 

-1 

-0.991 

0.9992 

1 


Table 3.2b 


Summary 

of the 

Comparison 

Results for 

Example 1.2 

s 


SS 

E b o’ b 0 

[ b i, b 2 ] 

[^2 



NG 

12 

16 

8 

IE --2 


Mad 

2E~2 

2E-2 

21E-3 



Mid 

QE-1 

IE -9 

0E~1 



NG 

12 

20 

12 

IE -4 


Mad 

22E-5 

29E-5 

16E-5 



Mid 

OE-1 

4E--6 

OE-1 



NG 

12 

24 

12 

1E~6 


Mad 

15E-5 

14E-5 

12E-5 



Mid 

IE -9 

IE -8 

IE-10 




Table 3.3a 


Choice of Stretching Points for Example 1.3 


8 

SP 

b o 

b l 

b 2 

b 3 

IE-2 

-1 

-0.98 

0.976 

1 

IE-4 

-1 

-0.992 

0.998 

1 

IE -6 

~1 

-0.9993 

0.999 

1 


Summary 

of the 

Table 

Comparison 

3.3b 

Results for Example 1.3 

e 


SS 

iv b i 

3 ■ 

[ b 1 ,b 2^ 

[b 2 



NG 

8 


16 

12 

IE-2 


Mad 

14E-3 


20E-3 

28E-3 



Mid 

OE-1 


4E-6 

OE-1 



IMG 

16 


16 

16 

IE -4 


Mad 

25E--5 


35E-5 

8E-4 



Mid 

OE-1 


2E-6 

OE-1 



NG 

16 


24 

16 

IE-6 


Mad 

33E-6 


19E-6 

9E-5 



Mid 

OE-1 


IE -7 

OE-1 





Table 3.4a 



Choice 

of 

Stretching 

Points for 

Example 

1.4 


£ 

SP 


b l 

b 2 

b 3 

b 4 

b 5 

IE 

-1 

-1 

-0.89 

-0.01 

0.01 

0.91 

1 

IE 

-2 

-1 

-0.849 

-0.007 

0.006 

0.87 

1 

IE 

-3 

-1 

CO 

« 

o 

i 

-0.006 

0.006 

0.84 

1 


Table 3.4b 

Summary of the Comparison Results for Example 1.4 


77 

SS 

[b 0> b i] 

Cb 1 ,b 2 ] 

[ b 2 9 b ^ 3 

[ b 3» b 4 ] 

[b 4 ,b 5 ] 


IMG 

8 

8 

8 

8 

8 

IE-1 

Mad 

24E-3 

* 

* 

* 

21E-2 


Mid 

IE-12 

* 

* 

* 

22E-7 


NG 

12 

8 

8 

8 

12 

IE-2 

Mad 

7E--4 

* 

JH. 

■& 

HE -4 


Mid 

OE-1 


* 

* 

12E-7 


IMG 

16 

8 

12 

8 

12 

IE-3 

Mad 

12E-5 

* 

-* 

* 

9E-5 


Mid 

OE-1 


J /r 

* 

OE-1 


■^indicates overflow in the numerical solution. 



Table 3.5a 


Choice of Stretching Points for Example 2 


e 

SP 

b 

0 

b l 


b 2 

b 3 

IE' 

-2 

0 

0.06 


0.95 

1 

1£' 

“4 

0 

0.005 


0.994 

1 

IE- 

"5 

0 

0.0002 


0.993 

1 


Table 3.5b 


Summary of 

the Comparison 

Results for 

Example 2 

-f-I- ss 

C b o ,b l-^ 

[ b !’ b 2 ] 

[bo ] 

NG 

8 

16 

8 

IE -2 Mad 

25E-3 

9E-2 

3E-3 

Mid 

0E~1 

2E -6 

OE-l 


NG 16 16 8 


1E~4 

Mad 

Mid 

28E-5 

OE-l 

32E-5 

3E~6 

4E-4 

OE-l 


NG 

12 

20 

8 

IE-5 

Macl 

31E-6 

25E-6 

2E-5 


Mid 

OE-l 

5E-7 

OE-l 





yu 

2 ~ 


1 



Figure 3.2 

Numerical solution to example 1.2 


—4 

,£=10 


-1 


0 


1 t 



Figure 3.3 

Numerical solution to example 


-4 

1 3, £=z\0 



3.4 Discussions 


Examples (1.2) - (1.3) exhibit boundary 
layers near both the end points t = + 1. The 
asymptotic solutions are constructed upto G( e^) 
accuracy, which can be accepted as fairly well 
approximations of the solutions to these examples. 
From the results in Tables 3.2 -3.3 , we see 

that the numerical solution compares well with the 
asymptotic solution. Again, as e decreases, the 
numerical solution compares better with the asymp- 
totic solution. 

Example 1.1 demonstrates the effectiveness 
of the hybrid method for interior layers. Here 
also as e decreases, the numerical approximation 
compares well with the asymptotic solution. 

In case of Example 1.4 it appears that the 
method does not work properly. But a more involved 
thought reveals the fact that the general solution 

t 2 t 2 

y(t) = (c-j^ + c 2 f e e dt) e E 

of the problem remains unbounded as e approaches 
zero. And in the numerical solution, we see that 
for smaller e's, overflow occurs; which is exactly 
what we expect. It is also observed that whenever, 



there is overflow in the numerical solution, the 
asymptotic solution is also too large. This 
demonstrates the efficiency of the numerical method. 
This also suggests that much of the qualitative 
natures of the solution might well be known by 
solving problems with the present numerical method. 

Example 2 illustrates the effectiveness of the 
hybrid method for higher dimensional problems. 

The figures indicate the absence of small 
oscillations in the numerical solution. In fact, 
the grid-generation helps to resolve the possible 
oscillations. It should also be noted that we have 
restricted the trapezoidal rule to the slow 

component which could otherwise produce the osc- 
illations. The numerical experimentation also 
tallies with the expected robustness of the method. 

In general, the thickness of the boundary 
layers or the interior layers are not known, 
analytically or asymptotically. But the choice of 
the stretching points (partition points) in the 
method indicates the numerical thickness of the 
layers. This provides us a necessary information, 
which could not be known otherwise. Unfortunately, 
this is not conclusive, since the grid-generation 



only chooses the stretching points by looking for 
the 'possibility' of the existence of the layers 
anc 1 not 'necessity'. However, after one gets the 
numerical solution, one can check, whether there 
is a rapid change of the solution in the concerned 
region; which would perhaps be conclusive for locat- 
ing the layers and determining their thickness(es ) 
numerically. 



Chapter 4 


SPECIAL MET H CDS 

4.1 Introduction 

As demonstrated in Chapter 3, the hybrid 
method solves linear systems of singularly perturbed 
two-point boundary-value problems accurately. The 
hybrid method is a general purpose method; and 
naturally it suffers from the demerits common to 
all the general purpose methods. The demerit is 
its lack of simplicity in computation. For simple 
homogeneous problems, without turning points and 
interior layers, the same procedure can be applied. 
Suppose that, from physical considerations, one 
knows that the problem at hand contains only the 
possible boundary layers, and for certainty there 
are no turning points or interior layers. So one 
naturally feels reluctant to proceed through the 
time-consuming grid-generation process and decoupl- 
ing. This is the beginning stage of looking for 
the special methods which would avoid such unneces- 
sary computation as well as save machine time. 

This chapter presents some special numerical 
techniques which avoid the cumbersome grid-generation 



and compute the numerical solutions directly. Note 
that the extra information the hybrid method was 
providing us is the layer-thickness. Even for a 
sub-class of problems such as homogeneous equations 
possibly with only boundary layers, for which we 
will apply the special methods, there is no way of 
determining the layer-thicknesses. This is rather 
a drawback of the special methods going to be pre- 
sented in this chapter. However, using asymptotic 
analysis and its procedure of solution, we can get 
a fairly good idea of this thickness, which we 
present in § 4.2. In $ 4.3 and § 4.4 two 
methods are suggested and also illustrated 

by means of some examples. The last section, 

§ 4.5 discusses the effectiveness and scope of 

applications of these special methods. 


4.2 Boundary Layer Thickness 

Consider the following singularly perturbed 
system of ordinary differential equations : 

x = A(t,e)x + B(t,e)y (2.1a) 

0 < t < 1 

ey = C(t,e)x + D(t,e)y (2.1b) 


with 



( 2 . 2 ) 


x i(°), x 2 (l), y 1 (0), y 2 (l) prescribed; 
where x = 

x x e m x , x 2 e m' 2 , Yl e ih" 3 , y 2 £iR l4 , and 

denotes d/dt, A, B, C, D are matrices of 
functions of t and e with appropriate dimensions; 
and e is a small positive parameter (0 < e << 1). 
Problems of this type have been well studied via 
asymptotic expansion techniques [41, 112]. Under 
the principal assumption that the eigenvalues of 
the (n^+n^) x (n^+n^) matrix D(t,0) have non- 
zero real parts, Harris [41 ] shows that such systems 
have an n ji +n 2 dimensional manifold of solutions 
which tend to the solutions of the reduced system 
(i.e., the system obtained from (l) by taking e = 0) 
as e -* 0. Moreover, if D(t,0) has k, 

0 < k < n 2 +n 4 stable eigenvalues, then there are k 
linearly independent solutions of (2.1) which decay 

U 4* / p 

to zero (like el , for some positive definite 
matrix ) as the stretched variable t = t/o -* oo, 
and there are n 3 +n 4 -k linearly independent solu- 
tions which decay to zero (like e” M 2 ( 1 " t )/ 6 for 
some positive definite matrix M 2 ) as the stretched 




Since the general solu- 


variable cr = ' (l~t)/e -» <*>. 
tion of (2.1) is a linear combination of these 
(asymptotic) solutions, any solution of (2.1) 
satisfies a three time-scale property. That is, 
such an asymptotic solution needs to be an additive 
sum of functions of t , t and a where functions 
depending upon t and cr represent the boundary- 
layer behaviours of the solutions at t = 0 or 
t = 1 and function dependent upon t satisfies 
the reduced system corresponding to (2.1). 

The question is exactly when one can neglect 
the functions of t and a to represent the solu- 
tion of (2.1) only by the reduced solution. In the 
literature, these regions beyond which the reduced 
solution represents the solution of (2.1) are 
termed as boundary layers and their measure, as the 
boundary layer thickness ( es ) . An inequality concern- 
ing this thicknesses was given in Hsiao and Jordan[46] 
for scalar equations of the type 

eu + p ( t ) u + q(t)u = f(t) 0 < t < 1 (2.3) 

with u( 0) , u(l) prescribed. 

In this section, an estimate of this thickness for 
systems of the type (2.1) is tried; and an inequality 



concerning the norm of exp(D( 0,G)th) is given, 
'th* being the thickness of the boundary layer(s). 

For our purpose, the following theorem will be 
sufficient [16,41,83]. 


Theorem 2. 1 

Assume that 


(H-l) 

( H-2) 
(H-3) 


(H~4) 


All eigenvalues of the matrix D(t,0) 
have non-zero real parts on [0,1J, 

D(t,0) is non-singular for t £ [0,1], 

The reduced system 

X = (A(t,0) - B(t,0)D” 1 (t,0)C(t,0))X (2.4a) 

X, 


X = 


X 


2~ J 


, X 1 (0) = x 1 (0), X 2 (l) = x 2 (l) 

(2.4b) 


has a unique solution. 

The initial boundary layer problem 


d„ 
d t 


" U 1 



I i i 

= D ( 0,0) 

n 

_ U 2_ 


L 2 J 


(2.5) 


has a decaying solution on 0 <. t < oo for- any 
initial value 0^(0); U 2 G B 4 ; 

(H-5) The terminal boundary layer problem 


6a 


V. 


V, 


= - D ( 1, 0) 


V. 


V, 


( 2 . 6 ) 



has a unique decaying solution on 0 < a < °° 

n o 

for any initial value V 9 (0); V 1 G3R , 
n 4 

v 2 e m 4 . 

Then, the system (2.1 ~ 2.2) has a 
unique solution of the form ; 

x(t,e) = X(t) + 0(e) (2.7a) 

y(t,e) = Y(t) + U(t) + V(cx) + 0(e) (2.7b) 

on 0 < t < 1, where the terms X, Y, U, V 
all have asymptotic power series expansions 
and the functions of t = t/e and 
a = (l-t)/e decay to zero as t , a - 

Remark 2.1 

Instead of assuming the existence-uniqueness 
of the solutions of (2.4 - 2.6) explicitly as in 
(H-2) -- (H-5). We could have taken some sufficient 
conditions to guarantee these [ 16 ] . 

The end-conditions to be supplied with (2.5) 
and (2.6) for getting a solution of (2.1) are 

U x (0) = y x (0) - Y x (0), U 2 (co) = 0 (2.8) 

and 


V 2 (0) = y 2 (l) - Y 2 (l), V x (») = 0. 


(2.9) 



But, since the solutions of boundary-layer problems 
(2.5) and (2.8), and (2.6) and (2.9) become insigni- 
ficant outsidide the boundary layer regions (due to 
their exponential decaying nature), it is appropri- 
ate to replace the conditions U 0 (“) = 0 and 

V 1 (°°) = 0 by U^T^) = 0 and V 2 .(^ 2 ^ = 0 for some 
x^ , o 0 < oo, respectively [109]. We require the 
solutions of the modified problems (taking x^ or 
cr 2 in place of <») to differ from that of the unmodi- 
fied ones only by 0(e) in magnitude (see equation 
(2.7)). We will use this requirement to have an 
idea of the thicknesses x^ and o 2 of the 
boundary layers near the end-points t = 0 and 
t = 1 respectively. For this purpose we change the 
notations slightly and give the modified and unmodi- 
fied solutions below. 


Let F = D(0,0) and U,V be the solutions of 
the unmodified and W, Z be the solutions of the 
modified boundary layer problems respectively. 


Let 




u x , w x e ih 1 

n o 

v x , z x e ir J 


n. 


U 2 ,W 2 e IR 2 • 


v 2 ,z 2 e ir 


n. 
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Then U, V, W, Z satisfy the following systems 
respectively. 


dU/dx = FU 

0 < T < oo 

(2.10a) 

U x (0) = Yl (0) - Y 1 (0) 


(2.10b) 

lim U 9 (x) = 0 

X-+ oo 


(2.10c) 

c’V/da = - FV 

0 < a < oo 

(2.11a) 

V 2 (0) = y 2 (l) - Y 2 (l) 


(2.11b) 

limV,(a) = 0 

a 


(2.11c) 

dW/dx = FW 

0 < X < T. 

(2.12a) 

W x (0) = y x (0) - Y 1 (0) 


(2.12b) 

W 2 (x 1 ) = 0 


(2.12c) 

dZ/d o = - FZ 

0 < a < a 2 

(2.13a) 

Z 2 (0) = y 2 (l) - Y 2 (l) 


(2.13b) 

Z,(a 9 ) = 0 


(2.13c) 


In the following, we try to find cut some sort of 
estimate of x^ by considering the systems (2.10) 
and (2.12). The results for o^ will follow similarly 
from the systems (2.11) and (2.13). 



A crude estimate can be obtained by taking the 
solution U of (2.10) and requiring that U 2 (e ) 
passes through an e -neighbourhood of W 2 (t) got 
from W(t), the solution of (2.12). But, since 
(2.10) is a two-point boundary-value problem, we 
are till now not in a position to solve it explicitly. 
However, we might use some sort of shooting to 
achieve this. So assume that 

(H-6) U 2 (0) = a < co f aem* 

such that lim U 0 (t) = 0 holds. 

Then, the solution of (2.10) is : 


U(t) 


Ft 


y x (o) - y 1 (o) 

a 


The requirement |U(x ^ ) | <_ e , for some 
norm | . j , gives us the following. 


Lemma 2.1 

Under the assumption (H~6) , t^ satisfies 
the inequality 


Ft. 

y x (o) - Y 1 (o)1 


e 1 

a 1 


1 

j 



F t 


~ E 1 (t) E 2 (t) 

E 3 (t) E 4 (t) 


(2.14) 


Moreover, if e 


(2.15) 



then a 


t 1 ^ U^CO E 3 (t)( Y (0) - Y,( 0) ) ] 

(2.16) 


Now, looking from a different angle we find 
that the difference f; = U ~ W satisfies a system 
of the same type as (2.10). The only difficulty 
being the end-conditions. So assume that 

( H-7 ) W 2 (0) = p < oo, p e m n4 

such that W 2 ( t 1 ) = 0(e) holds. 


Then £ 



n. 




m 


n y 




satisfies the following system : 


d^/dx = F £ 0 < t < t x (2.17a) 

(0) = 0, g 2 (0) = a - p (2.17b) 

The solution of (2.17) is 

t<* ) = * FT T ° " • 

_ a-p __ 

Or, using (2.15) we. can. write it as 

E 2 ( t) (a-p) 

E 4 ( x ) («-J3) 

As before, we set | ?( t^) I Js. ^ 2 » 6 obtain 


= 



Lemma 2.2 


Under the assumptions (H-6) - (H~7),t,j satis- 
fies 


j E 2^ t i^ ( a -P) 

I E 4 (' r 1 ) (a~p) 


(2.18) 


for some positive constant c 3 > ’where a is deter- 
mined by (2.16) and J3 satisfies : 


I E 3 (T 1 )(y 1 (0) -Y 1 (0)) + E 4 (x 1 )p | < c 3 .e (2.19) 

for some positive constant c 3 . 

It is of special significance to consider the 
degenerate cases arising from the trivialities 
n 3 = 0 and n 4 = 0. In these cases, the asymptotic 
solution of (2.1) - (2.2) satisfies a two-time scale 
property instead of the three-time-scales, that is 
because one of the boundary layers is lacking. 


Case 1 When n 3 = 0. 

Here, y = y 3 « Consequently the boundary condi- 
tions (2.2) are given by x^O), x^l) and 
The reduced system approximates the original one 
on [0,1). At t = 1, the boundary-layer correction 
V(cr) satisfies : 



dV/da = 

-FV 

0 < a < co 

(2.20a) 

V(O) = 

y( i) - 

Y(l) 

(2.20b) 

i im V ( a ) 

0-4.00 X ' 

= 0 


(2.20c) 

Setting 

|V(a 2 ) 

| <, c 4 . e , we obtain, 


1 Y Fa 2 

1 0 

(y(D - 

Y(l)) | < c 4 .e 

(2.21) 


for some positive constant c 4 .e 

Case 2 When n 4 = 0. 

In this case y = y^. Thus we only consider 
the boundary layer correction at t = O to obtain 
the estimate (similar to Case 1) : 

FT 

I e x (y(0) - Y( 0) ) | < c 5 . e (2.22) 

for some positive constant c^. 

Remark 2.2 

If both of n.,, n^ are equal to zero, then 
(2.1) is not a singularly perturbed system and we 
might not observe the boundary layer phenomena. 

For non-homogeneous equations (unlike (2.1)), the 
same procedure might be applied, but it would 
involve the terms which might need 



extra hypotheses to yield practically meaningful 
results . 

We finish this section by citing an example 
from O'Malley [83] where, the thickness is directly 
obtained . 

Example 1 

Consider the following TPBVP : 


y l = y 2 



with the boundary conditions, 

Yj_(0) = 0, y 2 ( 1) = (be 6 - pe^)/ ( be - \i e) 

where 6 = (-1 + Yl~4~ 6 )/2 , 

\i = (-1 - Y l-4 s )/2 . 

Writing the above system as a second order equation, 
we obtain, for x = y^ , x = y 2 » 

ex + x + x = 0 . 

The exact solution x(t,s) is given by : 
x ( t ,e ) = ( be - pe)” 1 (e 6t - e P-t ) . 

If we neglect 0(e) terms, we obtain, 



x(t,e) e - e , where e ~"t/ e represents 
the boundary --layer corrections. The thickness t^ 
of the boundary layer near t = 0 can directly be 
obtained by taking e "^l^ 2 <_ c^. for some constant 

c 6 a '^ our disposal. Choosing c^ = 1, we obtain 
t j_ 2. |e In e | . 


4.3 Boundary Value Technique 

In the introduction to this chapter ( § 4.1), 
we have remarked that some special numerical methods 
would be developed for systems of two-point boundary- 
value problems where there occurs no turning points. 
In the last section ( <§ 4.2) we found out the thick- 

ness of the possible boundary layers (in fact, the 
exponential of the thickness). This helps us to 
choose the points t^, t 2 in [0,lj properly so that 
the boundary layers might be expected to lie only 
within [0,t 1 ] and [t 2 ,l]. In this section, we pre- 
sent such a method which will be discussed only for 
linear homogeneous TPBVP's. The non-homogeneous 
case can be handled similarly. To this end, we con- 
sider the following problem (see equation (2.1)) 


x 


sy 


Ax + By 
Cx + Dy 


0 < t < 1 


(3.1a) 



with boundary conditions 


x 1 (0) = a , x 2 (l) = b , Yl (0) = c, y 2 (l) = d (3.1b) 



A(t,e) , B(t,e), C(t,e), D(t, e) are 

smooth matrix functions with appropriate dimensions 
and a, b, c, d are constant vectors, denotes 

d/dt and e is a small positive parameter. 

Though , .this class of problems is a 

very restricted one, this is worth paying an atten- 
tion due to its appearance in many physical applica- 
tions (e.g. see Chapter 6). We assume that there 
is a unique bounded solution to (3.1) and that 
there are no turning points or corner layers or 
interior layers in the solution to (3.1). 


These assumptions lead, us to a simplified grid- 
generation process in the hybrid method. In fact, 
if we could determine t^ and t 2 , the only pos- 
sible stretching points, (i.e . , where a finer grid is 
required), we need only some grid with a modest 
grid length so that the numerical methods we might 



like to use are accurate . Note that we also need 
to determine the solution values at t^ and 
In that case, we divide the interval of integration 
[0,1] into three sub-intervals [0,t 1 ], [t x ,t 2 ] and 
[tgjl] and solve the stretched equations (according 
to the stretching parameters) on these sub-intervals 
independently. In the following, we describe this 
technique in detail and in a step-wise manner. 

S tep 1 

Select points t^ and t 2 in agreement with the 
inequalities of the form (2.18). Note that, usually 
the thickness of the boundary layers are 0(e| In e|). 
Determine the values of x and y at t^ and t^ using 
asymptotic expansion. To this end, let 

CO CO 

x = S x / • \ ( t ) 1 , y = Z y m (t) e 1 (3.2) 

i=o ^ ; i=o v ' 

where the subscript i within paranthesis shows the 

formal terms x^^(t)'s as t-dependent co-efficients 

of e 1 and etc. on substituting (3.2) in (3.1), the 

resulting equations can then be readily solved as 

initial value problems for x^^(t) and y^)(t) *, 

using analytical, or numerical means [58,28,112,82]. 

Once we get the values of x(t^), y(t.) for j = 1,2, 

r\ 

either with 0(e) or 0(e) expansions, we switch 



over to the next step, wherein we obtain a solution 
to the TPBVP (3.1). 

Step 2 

We scale the boundary layers by the factor e . 
So denote 

T(t) = t/e , cr( t) = ( 1-t )/e (3.3a) 

T 1 =T ( t i)» (3.3b) 

Since, the boundary layers are guided by the t ~ 
and o-scaled equations, we resolve (3.1) into the 


following TPBVP' s s 
Left boundary layer : 

dx/d t = e(Ax 4- By) (3.4a) 

dy/dr = Cx + Dy (3.4b) 

on 0 < t < with boundary conditions : 

Xj(O) = a, Yi(O) = c (3.4c) 

x 2 (t 1 ) = x 2 (t 1 ), y 2 ('C 1 ) = Y2 (t l ) (3.46) 


Note that in (3.4d) above, the right hand side 
terms x 2 (t 1 ), y 2 (t 1 ) are as found out in ste P 



Original system 

The equations (3.1a) with t £ [t^jt^]; the 
boundary conditions being 

x l(ti) , Y]_(tj_) and as obtained 


in step 1 . 

Right boundary layer : 

dx/do = - £ (Ax + By) (3.5a) 

dy/da = - Cx - Dy (3.5b) 

on 0 < o < Op , with boundary conditions : 
x 2 (0) = b , y 2 (0) = d (3.5c) 

x l(a 2 ) = X x (t 2 ) » y x (cr 2 ) = Yj_(t 2 ) (3.5d) 


Then the above three systems are solved numeri- 
cally and are combined to give a uniformly valid 
solution on [0,l], The left boundary layer and 
the right boundary layer are solved by using Euler's 
Implicit rule. We apply this rule, since the left 
and right boundary layers are stable in t and -t 
respectively 5 the reverse stability of the right 
boundary layer is taken care of in formulating the 
system ( 3 . 5 ) by the transformation a = (l~t)/e . 

In the middle of the interval of integration, i.e. , 



the original system as written above is solved by 
using trapezoidal rule. This agrees with the 
hybrid method, since on this interval the solution 
is slowly varying. 

Remark 3.1 

The boundary conditions in (3.4 - 3.5) and for 
the original system are chosen somewhat arbitrarily. 
This way of choosing the components x^, x 2 , etc. 
and recalling their values at t^, t 2 comes in 
accord with the asymptotic choice (see the cancella- 
tion rules in O'Malley [82])* For example with 
(3.4a -- 3.4b) we have to supply conditions at t = 0 
and t = First, at t = 0, in (3.1b) we have 

the values x-^(O) and y^(0). They are retained. 
Then we need the values of x 2 and y 9 at the 
points 0 or, t^. But they are obtained only at 
in step 1. Hence the conditions (3.4c - 3.4d). 
Analogously boundary conditions in (3,5) are chosen. 
Then, the remaining boundary conditions are supplied 
to the original system. 

Remark 3 .2 

Asymptotically, it is instructive to solve the 
reduced problem on the interval [t^,t 2 J. But this 



creates further problem as to which boundary condi- 
tions are to be chosen; since the dimension of the 
reduced problem is n 1 +n 2 compared to ' (3.1) 
which has the dimension n 1 +n 2 +n 3 +n.. The other 
reason is that, a better approximation is expected 
if we solve the original system instead of the 
reduced system. 


Step 3 


In choosing t-^ and t 2 , note that we could 
only determine them upto an 0(e | In e | ) . So it is 
necessary to check (numerically) whether the choices 
are judicious. So let X(t) = j^y^j » let X(t^) 
and X(t^) denote the approximations (as obtained 
from the numerical solution in step 2) of X(t) 
computed from equations (3.4) and the original sys- 
tem; let X(t^)_ and X(t 2 ) + denote that of X(t) 
at t = t 2 obtained from the original system and 
equations (3.5) respectively. Then use the follow- 
ing criteria. If, 


I X(t.), » X(t.) | < 6 , i = 1,2 (3.6) 

X ~r X *— 

is satisfied where jaj denotes the maximum of the 
absolute values of the components of a and 6 is 
a prescribed tolerance limit (of 0(e) or 0(e )), 
then we stop; else, we choose different t^, ^2 



still closer to t = 0 and t = 1 and proceed 
again from step 1 to step 3. In fact, we solve 
first by taking t^ = 10 e j In e j and etc. and 
compare the norms in (3.6) for different iterations 
by making t^ smaller (similarly for t^), so that 
we can restrict ourselves to find better choices of 

t l» t 2* 

Remark 3.3 

Since, the technique above heavily depends 
upon the choice of t^ and ± 2 ) and since the 
values of x and y at these points are determined 
using the asymptotic expansions in the boundary 
layers; we should* not increase but decrease t^ or 
l-t 2 in subsequent iterations. The name 'Boundary 
Value Technique' is given for this reason. 

From among several test examples, the following 
is selected to illustrate the technique. The boun- 
dary conditions, are also specified for an even 
illustration. 

Example 1 

Consider the following second order problem 
(cf. [58]) 

ey + ( l+2t )y + 2y = 0 0 < t < 1 (3.7a) 

y (0) = 0, y(l) = 1. 

Writing this as a first order system, 


we obtain : 



similarly the boundary conditions are supplied. The 
similar procedure is followed for g = io~" 4 , 10”^ 

and 10~ 6 . 

The results for e = 10“ 3 , 10" 4 , 10~ 5 , 10~ 6 as 
obtained in the third iteration are summarised in 
Tables 4.1. In Table 4.1b, we give the comparison 
results , obtained by comparing the numerical solu- 
tion produced by the boundary value technique with 
the asymptotic solution 

y(t^ ) =3 ( -j~ - e“ e ) +e[ + e“ i ) 

t 2 t 

+ 6( — * - e" e ) - 3— e _ e ] + 0(e 2 )(3.9) 

(l+2tr e 2 

In Figure 4.1 the solution y(t) is shown when 
-5 

e = 10 . Additional headings in the tables are as 

follows : 


e_ = 

max | 
t6BL 1 

U A (t) 

- u N (t)| 

(3.10a) 

II 

o 

Q) 

max | 
t£BB 1 

u A (t) 

I 

c 

:2 

<- 1" 

(3.10b) 

e + = 

max | 
t£3R 1 

u A (t) 

- u N (t)| 

(3.10c) 

where 

BL, BR 

and 

BB are 

the left boundary layer. 


the right- boundary layer and the region beyond the 

boundary layers respectively. u A (t) and u N (t) 

are the asymptotic and the numerical solutions respec- 


tively. 



Table 4.1a 


Choice of the Points t 1 and t 2 for Example - 1 


e 

CO 

f 

O 
• — ! 

1C' 4 

10“ 5 

10- 6 

t i 

4xl0” 3 

3.2xl0“ 4 

4.4xl0~ 5 

5xlO“ 6 

i 

CM ! 
•p 1 

no choice 

in 3rd 

iteration 



Summary of the 

for Example 1. 

Table 4.1b 

Comparison Results (3rd iteration) 

6 10~ 3 

10“ 4 

10" 5 

10~ 6 

e_ 0.2426E-1 

0. 6721E-2 

0.1492E-3 

0.3720E-4 

e Q 0. 3127E-2 

0. 7823E-3 

0. 2129E-3 

0.4521H-4 


e + 


Ho right boundary layer 


0 


1 


t 


Numerical 


Figure A A 
solution to 


-5 

example 1 , E =*=10 


4.4 Cutting Point Technique 

In § 4.3, we have presented the Boundary 
Value Technique to solve linear homogeneous singu- 
larly perturbed systems of two— point boundary value 
problems. The choice of the points t^, that 

separate the inner regions and the outer regions 
(as they are usually called) is made partially by 
using the estimates in § 4.2 and modifying it appro- 
priately by the successive iterations. Note that 
the estimates in $ 4.2 yield an inequality of the 
form 

c je In s | t^ for some constant c. 

First we choose greater values of t^ and reduce 
it to such a value which would separate the inner 
regions and the outer region accurately. 

In this section, we will follow the reverse 
trend, we first start from a value of t^ , (say, 
equal to e) and increase it till we get the right 
value for the thickness of the boundary layer. Con- 
sequently during the computation the inner regions 
are made bigger and bigger in subsequent iterations 
so that the points t-^ and tg 9° inside the outer 
regions where asymptotically, the reduced solution 
is a 'good* approximation to the exact one. 



For this reason, the first step in the Boundary 
Value Technique is replaced by the following which 
we write as step 1 (for the present method, the Cut- 
ting Point Technique). 

Consider the equations (3.1). The Cutting Point 
Technique is described in a step-wise manner below. 

Step 1 

Choose [0»l] close to 0 and 1 res- 

pectively i.e., in the neighbourhoods of the points 
0 and 1 of radii of 0(e|ln e|). For the evalua- 
tion of x(t i ), y(t^), i = 1,2* take the reduced 
system corresponding to (3.1). The reduced equations 
are ; 

x = Ax + By (4.1a) 

0 = Cx + Dy (4.1b) 

with the boundary conditions : 

x x ( 0 ) = a , X 2 (l) = b. (4.1c) 

Here, we assume that (4.1) has a unique bounded 
solution. Then evaluate x(t^) and y(t^), i = 1,2 
from (4.1) either by analytical or by numerical 


means . 



Step 2 


Step 2 of the Boundary Value Technique serves 
as step 2 of txhe Cutting Point Technique also. 

Step 3 

We require that (see step 3 of Boundary Value 
Technique in § 4.3) 

I X + (t i ) ~ X_ (t.) | < 6 , i = 1, 2 (4.2) 

If this requirement is not fulfilled, we choose 
and tp farther from 0 and 1, i.e., we increase 
the thicknesses t^-0 and l-tp and go to step 1 
i.e., we choose t-^ greater than and t^ less than 
that previously chosen and iterate step 1 to 3. 

Among several test examples, the following 
example is chosen to demonstrate the method. It is 
noteworthy that, the method when applied to 
Example 1 in § 4.3 yielded an approximation in the 
4th iteration which matches with the 3rd iteration 
of the Boundary Value Technique upto four significant 
decimal places. 



E xample 2 

Consider the second order constant 
equation 

ey + y + y = 0 0 < t < 1 

y(O) = 0 , y(l) = l/e . 

Or, in system form, 

* 

y = z 

m 

ez = -y - z 

on 0 < t < 1 with boundary conditions 

y(O) = 0 , z(l) = l/e . 

The exact solution of (4.3) is : 

— \i t — Xot 

y = c^e l + c 2 © ^ 

with ^ = (1 - yl-4ej/2 e 

■h 2 = (l + YT-T4e)/2 e 

C 1 = (e C 3 )" 1 , c 2 = -(e c 3 ) 
c 3 = ^ ^ e- X l . 


coefficient 

(4.3a) 

(4.3b) 

(4.4a) 

(4.4b) 

(4.4c) 

(4.5) 


The results are 


summarised in Table 4.2 and the 



solution is shown in Figure 4.2 for e = 1CT 4 . 

Below, the tables are given for the 3rd iteration, 
where there is no need to choose as it appears 

from the first iteration that there is no right 
boundary layer. 



Table 4.2a 


Choice of the Point t-^ for Example 2 


e 

co 

f 

O 

< — i 

10" 4 

10" 5 

10~o 

i — i 

-P 

3xl0“ 3 

4xl0" 4 

9xl0~ 4 

5xl0~ 6 



Table 4. 

2b 


Summary of the Comparison Results ( 3rd 

iteration) 

for Example 2. 




e 

CO 

l 

o 

10" 4 

10" 5 

10" 6 


0.4321E-2 

0.2372E-3 

0.4998E-4 

0. 1217E-4 

e o 

0.5789E-2 

0.2319E-2 

0. 1986E-3 

0.4217E-4 


Table 4.2a 


Choice of the 

Foint t^ 

for Example 2 


e 

10" 3 

1CT 4 

10~^” 

10~ 6 

h 

3xlO" 3 

4xlCf 4 

9xlG" 4 

bxlCf 6 


Table 4.2b 

Summary of the Comparison Results ( 3rd iteration) 

for Example 2. 



0.4321E-2 


e o 


0.b789E-2 


0.2372E--3 0.4998E-4 0.1217E-4 
0.2319E-2 0. 1986E-3 0.4217E-4 



Figure 4-Z 

Numerical solution to example!, £-10 


4.5 Discussions 


For the results of the numerical experimenta- 
tions with the special methods, namely, the Boundary 
Value Technique and the Cutting Point Technique, we 
have seen that there is practically, a very little 
difference in the solutions obtained by the computa- 
tional and asymptotic methods. This shows the reli- 
ability of these special methods. One should also 
look at the sub-class of problems these special 
methods tackle. These methods, we hope, might be 
applicable to the non-homogeneous linear systems 
but without interior layers and etc. We shall try 
latter to apply these methods for non-linear systems 
(see chapter 5). However, for non-homogeneous cases, 
the estimates in .§4.2 are not sharp. So when we 
try to fix the size of the boundary layer regions 
(by trial and error) via the iterations, we must 
not go in one direction, i.e., we must not only 
decrease the boundary layer thicknesses throughout 
the computation nor that we must only increase it 
throughout the computation. That is, our choice 
of the thicknesses would oscillate from iteration 
to iteration. This way of trying for the thicknesses 
(determining appropriate ^ and t 2 ) might be 



incorporated in homogeneous cases also; which, 
of course, will demand more machine time. 

Another consideration in this direction is to 
use other methods, like Simpson's rule (e.g., Reddy 
and Kadalbajoo [51], Kellogg and Tsan [57]) in place 
of Euler’s Implicit, and Trapezoidal rules. We 
have used these due to the analysis of the hybrid 
method in chapter 2 . 

Specifically for homogeneous systems a new 
technique would result, if we can eliminate some 
of the concerned variables to reduce the system 
into a second order system of lower dimension. 

This might be done in case the dimension of the 
system is (2m+2n). To be precise, consider the 
following system (which often arises in Optimal 
Control Theory, see Chapter 6). 

x = a^x + + + B^q 

ey = A 3 x + A 4 y + B 3 p + B 4 q 

p = C ] _x + C 2 y + D 1 p + D 2 q 

eq = C 3 x + C 4 y + D 3 p + D 4 q 

on 0 < t < 1 where x,p£m m , y.qem" and all the 

A.(t f O, B,(t, a), qtt, e) and 
1 X 


(5.1a) 
(5.1b) 
(5.1c) 
(5. Id) 


other matrices 



D^(o,e) are of appropriate dimensions. For elegance, 

let 




and E = 

I- being the identity matrix of order j. Then, 

3 

the system (5.1) might be rewritten by 




. . , , + q io non-sinaular throughout [0,l], 

Assuming that B is non smyu 



from (5.2a) we obtain 


P 

q 



x 

y 



(5.3) 


with 





9 


Ff — (B^ - ti^B^ B^) , 

F 2 » b“ 1 b 2 f 4 , 

Putting (5.3) in (5.2b) 
system 


F 4 = (B 4 - B 3 B' 1 B 2 r 1 , 

F 3 = b“ 1 b 3 f 1 . 

, we Obtain the second order 




dt‘ 


,.y. 


® _lA dt 


Ly. 


+ (A - BB~ i A) 


x 

. Y. 


(5, 


The same form could be obtained if instead we have 
D a non-singular matrix, where we can eliminate 


X 

y 


instead of 


P 

q 


tions can be obtained for 


. The boundary condi- 
by using (5.3). 


x 

y 


Klomq one might find which ever 
For particular problems one m y 



way is convenient. The advantage in doing this is 
that one can apply the results for second order 
systems (see, e.g., Howes [45], Eckhaus [28] and 
Engstrom [29]). Also, the interval of integration 
can be taken as any bounded closed interval [t ,t £ ] 
in place of [0,l]. 

In summary, we are justified in trying to find 
a solution in the left boundary layer with the 
stretching transformation x = t/ e and in the right 
boundary layer, with o = (l-t)/e , since the boun- 
dary layers are stable in time t and reverse time 
-t respectively. Again, the approximation in the 
outer regions (i.e., beyond the boundary layers) by 
the proposed methods are expected to be better than 
the asymptotic solutions, since in the proposed 
methods, the original system is solved, whereas in 
the asymptotic techniques, the solution in the outer 
region is represented totally by the reduced solution. 
One is justified to approximate the solution, in this 
outer region by the reduced solution, when e is 
small enough, i.e., an infinitesimal (see e.g., [73] 
and etc.* also a beautiful counter example in 

Kokotovic and Yackel [61], ™ he n 8 ’ °- 1 ^ • But in 
practice, e is not an infinitesimal; whence it 



amounts to make the proposed methods preferable. 

Also, since the thickness of the boundary layers 
is not yet determined with analytical or asymptotic 
means for the systems in hand, and asymptotic methods 
available to solve this class of problems almost 
avoid this question, it is suggestive that we try 
to get a solution by these special methods, where 
we would be able to know more accurately about the 
thickness of the boundary layers (rather than just 
0(e ) or 0(e | In e | ) ) . To note a few, among the merits, 
the methods do not require much skill necessary in 
the conventional methods (e.g., matching the coeffi- 
cients in the asymptotic expansions), it does not 
require very fine mesh size; the methods require 
very less problem preparation for computer adapta- 
tion , these are also efficient and accurate; the 
numerical experiments demonstrate this fact. It 
might also be observed that an accuracy predicted 
could always be achieved with comparatively less 
computational effort. 



Chapter 5 


NONLINEAR SINGULARLY PERTURBED TWO-POINT 
BOUNDARY -VALUE SYSTEMS 

5.1 Introduction 

So far we have presented some efficient and 
accurate methods for linear systems of singularly 
perturbed two-point boundary-value problems. In 
this chapter, we consider the corresponding non- 
linear problems of the form : 

x = f(x,y,t,e) (1.1a) 

0 < t < 1 

ey = g(x,y,t,e) (lilb) 

with boundary conditions 

B(x(0), x(l), y(0), y(l),e ) = 0 » (1.1c) 

x £ ]R m , y eiR n , f, 9 are smooth nonlinear func- 
tions; (1.1c) represents the boundary condition, 
which when linearized would provide m+n linearly 
independent conditions. We also assume the unique- 
ness and boundedness of the solution to (l.l) on 
the interval [0,1]. For numerical solution of 
(1.1) one often thinks of applying usual methods 
([55-56]). The experience with the linear case shows 



that the usual two-point boundary-value methods 
are not expected to work satisfactorily. Whatever 
the methods are, one should try to find out esti- 
mates of the type discussed in Chapter 2. But 
this certainly aims at too much of a difficult task, 
as the linearity, the most important of all, is 
lacking. 

This leads us to consider the other approach, 
i.e., using iterations. The procedure is that we 
linearize the problem (1.1) first and then find out 
numerical solution by the methods developed in the 
preceding chapters. Though this is quite plausible, 
the main disadvantage of this approach is in terms 
of cost. It approximates the system (l.l) by a 
sequence of linear systems* all of which then must 
be solved numerically. However, experience with 
some specific problems of the class defined by (l.l) 
shows that we might follow this approach. 

In § 5.2 such an iteration process, i.e., a 
modification of the familiar Newton's iteration is 
proposed. Then a model example is considered in 
§ 5.3, and the method thus developed is applied 

to this example. And in the last section, the 
remarks on the efficiency of the methods are made. 



5.2 Modified Newton's Iteration 


For notational convenience, we consider the 
nonlinear autonomous equations corresponding to 
Equations (l.l); we also supress the explicit 
appearance of e in (l.l), for the time being. 

So we consider the system 
* 

2 = F(z) 0 < t < 1 (2.1a) 

with non-linear boundary condition whose lineariza- 
tion yields 

B z(O) + B 1 z(l) = B 2 (2.1b) 

representing k linearly independent boundary 

conditions* k the dimension of z. 

* 

The Newton's iteration is defined in Theorem 
(6.6) of Chapter 1 where also the quadratic conver- 
gence is stated. However , for practical purposes, 
the Frechet derivatives or the G-derivatives do not 
play very important role. So a modification of the 
iteration can be carried out by replacing the Frechet 
or G~derivatives by the corresponding Jacobian dF. 

To this end, we use the generalized Taylor’s theorem 
for operators (Theorem 6.3 of Chapter l) to obtain 



I |P(z) - PU 0 ) - dP(z 0 )(z-z 0 )|| 



by taking n - 2 in Theorem (6.3) of Chapter 1, 

Note that dP(z ) and d 2 P(w) denote 6P(z )/dz 

u 0 0 

and d P(w)/dw' respectively. Newton's iteration 
might be written (for this case) by : 

z n+l = z n + z n n = 0, 1, 2, (2.3a) 

where z n satisfies : 

d p (z n ) z n = - P(z n ) 5 (2.3b) 

The subscripts n in (2.3) denote the nth 
approximation of z, the solution of the equation 

P ( z ) = 0 (2.4) 

It follows that the iteration (2.3) converges quadra- 
tically near z; since, (2.2 - 2.4) implies that 

I U n+1 - z 1 | < I j [dP(z n )3 1 | I . 

sup , \ | 1 d 2 P(w) j |.j |z n -z|J 2 , 
weL(z n » z ) 

assuming that P(w), [SP(w)]- 1 and S 2 P(w) are 



bounded in some neighbourhood N(z,r) with P(z) = 0. 
Now, our aim is to apply the iteration (2.3) to (2.1). 

Here, instead of applying the iterations above 
directly by letting P(w) ~ w - F(w), we will go 
a little bit further. Suppose that the nth approxima- 
tion z n is already determined. We want to find out 
the (n+l)th approximation z n+1 « 8u-fc z n+l satisfies 
(2.1a) which means : 


z n+l 


F < z n + 1> 


(2.5) 


Our scheme is to approximate F( z n+ i) by the method 
of operators, i.e. , we let 

FU n+1 ) ~ F(z n ) + K(z n )(z n+1 - z n ). (2.6 


Note that, formally we are not justified in 
approximating this way. Assuming that this linear 
ization is acceptable, we obtain : 


W ■ F(z n> 

Subtracting z 


z n+l ” Z n 

Denoting z n = 


+ dF(z n )(z n+1 - z n ) 

from each side of (2.7), 
n 

dF(z n )(z n+1 - z„) +(F(z n ) 

. 7 _ z , we obtain 

■ z n+l n ’ 

z„ + (F(z n ) - Z n ) 


(2.7) 


we get 
z n ^ ( 2 * 8 ) 


(2.9) 



Now, (2. 3a) and (2.9) define the iteration process, 

we are interested in . 

Rem ark 2.1 

'•* e have avoided the second derivatives that 
might have come up in considering the operator 
P(w) = w - F(w) directly. The second point is 
that, one should not neglect the term F(z ) - z n 
in (2.9), no matter how much temptation it might 
offer, because, z n might not be approximating the 
exact solution. 

In the following we shall show that the conver- 
gence of the modified Newton’s iteration (2.3a) and 
( 2 . 9 ) is also quadratic. We begin with an initial 
value problem (2.1a - 2.1b) and 

z ( 0 ) = a < 00 . (2.10) 

Very often in Mathematics, it happens that the 
same result is proved again and again. This seems 
futile but it is not. Our proof of convergence 
below begs the same excuse. 

Lemma 2.1 

Let w(t)€B^ be the solution of the initial 


value problem 



w = 5F(y)(w~y) + F(y) 


0 < t < 1 (2.11a) 


w(0) y (0) (2.11b) 

y(t) being a known function of t. Assume that 

(a) x£ lTv,w) II^WII ■ 0(1) 

(b) exp( / I |6F(y(5)) 1 |d£ ) = 0(1) , 0<_ rj < t < 1 

T) 

Then, | | F ( w ) - w| | < = 0 (g<t<l I I ^ ( y ) — y | j ). 

P roof (2.11) implies ( [94] , Sec. 33) that 
j |w(t) - y(t)| | 

< / exp( / 1 |dF(y(g )) | |n§) | |F(y(m)) - y(n)||dTj 

O T) 

< sup j |F(y) - y| | / exp( / I |dF(y(|))| tdDdr) 

0<t<l 0 ^ 

< sup | |F(y) - y j I sup (t exp(}||S F ( y (6))||dS)) 

~ o<t<i ri 

X * 0 # f 

| |w(t) - y(t) 1 1 < K 2 = 0(^P< 1 l|F(y) -yli). 



Again, 

| |F(w) « w | j 

= j|f(w) - f (y ) - dF(y)(w-y)|| 

«E U lfy,«) H * 2 F(x) H 

from which the result follows. 

Lemma 2.2 

If z(t)£ B^" is a solution of 
z = F(z) 0 < t < 1, z (0) given ; (2.12) 

and w(t) £ B is a solution of 

w = dF(y)(w-y) + F(y) 0 < t <, 1 (2.13a) 

w(O) = z(O) (2.13b) 

with a given function y in a neighbourhood 
N(z,r), r > 0 of z and conditions (a) -- (b) in 
Lemma 2.1 are satisfied, then 

| | z--w | | < K = U (|iz-y|| 2 ). 

Proof 

Let s = z-w 

and G(y,z) = F(z) - F(y) - dF(y)(z~y) 

6.3 of Chapter 1, we obtain 


then from Theorem 



Also, from (2.12) - (2.13) we obtain 
s = dF(y)s + G(y,z), s(0) = 0. 

Then, as in Lemma 2.1 

ll s (t)|| < / exp( / | |dF(y(g))| |d | ). 

o r\ 

• I |G(y(tj)) z( r)> | | dr) 

< 2 x S eL(y,z) 1 1 d2p ( x ) | |t exp(/| JdF(y( ; g) ) j |dg) 

* 1|z(ti) " y(Ti) i * 2 

which proves the lemma. 

Immediately, we get the following theorem which 
states that the solution of (2.1) with an initial 
value z ( C) obtained by iteration has quadratic 
convergence. 

Theorem 2.1 

Assume that the conditions (a) — (b) in Lemma 2.1 
are satisfied, then the solution z*(= z n ^ 
obtained by the following iteration scheme : 

z n = dF ^n } \ + (F(z n } " z n } 


(2.14a) 



0 


(2.14b) 


z 


n 


■'n+1 


= z n + z 


n 


(2.14c) 


suitably chosen with z q (0) = z(O) (2.14d.) 


converges guacSratically to the solution z(t) of 
the equations 

m 

z( t) ~ F(z) 0 < t < 1 (2.15a) 

z ( 0) given. (2.15b) 

• 

Moreover, if ||F(z q ) - z Q J | < 1, then (2.14) 

solves (2.15) with | J F ( z n ) - z n |j converging 
quadratically to zero. 


Proof 

In Lemma 2.2 put w = z n+1 » y = z n to get 
the algorithm (2. 14). We also obtain 

liz-z n+1 || < K = 0 ( | |z - z n | | 2 ), proving 

the first part of the theorem. By taking w = z n+1 , 
y = in Lemma 2.1, we obtain 

| |F(z n+1 ) - z„| < Ki = 0 (||F(z n ) - z n || 2 ) 

which proves the second part. 

Instead of the given initial value z(0) 

that a separated set of boundary conditions 


suppose 



is given, e.g., 


Zj_(0) = p, z 2 (l) = Y> 2 



(2.16) 


xhen, also the lemmas hold, since in that case, we 

t 

replace the integrals of the type f <p ( . . . ) by 

o 


t 



f %(...) 

0 

t 

where cp = 

~ *1 " 

/ V (...) 

_ 1 


„„ 


However, for the TPBVP (2.1) we assume that it 
possesses a unique bounded solution. Then, since 
z(t) is bounded, we have, in particular, z(O) = a < ». 
Then the Theorem 2.1 will hold. We note this down 
in the following.. 


Theorem 2.2 

If the solution of (2.1) is unique and bounded 
and F satisfies the conditions (a) - (b) in Lemma 2.1, 
then the iteration scheme (2.14a - 2.14c) with z q 
satisfying (2.16) solves (2.1a~c) with quadratic 
convergence . 

Remark 2.2 

The condition z n = 0 in (2.14) might be 



replaced by B q z n (o) + B^l) = c where 
B q z( 0) + B 1 z(l) = B 2 is to be satisfied by the 
initial approximation z Q . This way of replacement 
should be carried out to solve (2.1) numerically. 

Note that we have considered only autonomous 
systems. However, it is standard to extend the 
theorems also to non-autonomous equations. Suppose 
that instead of (2.1) the following equation is con- 
sidered 

z = F ( t , z) . (2.17) 


Let w 


t 

z 


G(w) = 


1 

F(t,z) 


Then equation (2.17) is transformed to 


w 


G(w) , 


which is autonomous. Similarly we add t(Q) = 0 

with the boundary conditions and then apply Theorem 2.2. 

To solve the linear systems in the iteration 
steps (in (2.14)) we apply the hybrid method devel- 
oped in Chapter 2. The other methods, i.e., the 
special methods of Chapter 4 might also be applied* 
but in that case, one must see that there are no 
interior layers in the original non-linear problem. 



5.3 A Model Nonlinear Example 


In the previous section, we have discussed 
a modification of the Newton's iteration we intend 
to apply for the nonlinear singularly perturbed 
two-point boundary -value problems. The systems we 
are concerned with are of the form : 

x = f ( x , y , t , e ) (3.1a) 

« 

ey = g(x,y,t, e) (3.1b) 

on 0 < t < 1 with boundary conditions 
B(x(0), x ( 1 ) , y(0), y(l), e) = 0 (3.1c) 

with the usual conventions. 

Writing z = 

(3.1a - 3.1b) can be reduced to 

z =-- G(z ) (3,3 

Our aim is to write the linearization directly for 
(3.1). For this purpose, we first compute 5G. 

df/dw df/dy 

e - 1 dg/5w s" ld g/ d y^ 



x 

t 

y , 


G = 


f 

1 


- 1 , 


(3.2) 



where w = 



let w n = w n+1 - w n and y n = y n+1 - y n . Then 

the linearization (2.14) reduces to : 

For n = 0,1,2,. . . 

* 

”n = ^("n-Yn^n + 37 (w n> y n )y n + (f < w n’ y n ) " V 

(3.4a) 

* 

s y n = t^V^n + t^'V y n )y n + ( 9 (w n’ y n ) - y n } 

(3.4b) 


w, 


n+1 


= w n + w. 


n 


y n+l = y n + y n 


(3.4c) 

(3.4d) 


w and y are initialized so as to satisfy a 
o ' o 

linearization of (3.1c), i.e., 



1 

0 

o 

X 

1 


-x 0 (l)- 

B 

0 

r 

■< 

o 

o 

5 

+b i 

y„U) 


(3.4e) 


with t Q (0) = 0. 


( 3 . 4f ) 


Note that here we might also choose the initial 
functions so that 



i I F ( z o ) " z 0 I I ^ obtain the results of 

Lemma 2.2. Now, (3.4a - 3.4b) and (3.4e - 3.4f) 
are the linear systems we handled in Chapter 3. 

So we apply the hybrid method to solve the systems 
(3.4a — 3,4b) and (3,4e - 3.4f) for each n. As 
experience shows it only takes 3 to 4 iterations 
to obtain an accurate numerical solution. However, 
the choice of the initial approximation z Q is 
significant in reducing the number of iterations 
for achieving certain predicted accuracy. 


The purpose of this section is to apply this 
numerical technique to a model non-linear example 
taken from Kevorkian and Cole [58], where almost 
all the varied nature of the solution are present. 


For computational point of view, we add another 
step to the above iteration and solution by the 
hybrid method. Note that, it is not yet clear, 
where to stop the iteration. To this end, we 
require an inequality of the type 


to hold, where 6 is a tolerence limit, whicn we 
prefer to take CXh*"), where h is the maximum 
grid length (matching with the quadratic convergence 



of the Newton's iteration and with the accuracy of 
the hybrid method). 

The example we want to study is : 

ex + xx - x = 0 0 < t < 1 (3.6a) 

x(O) = A, x ( 1 ) = B . (3.6b) 

x(t) £ IR ; A and B will be specified later. 

Equation (3.6a) is invariant under the transforma- 
tion : 

x - — > -x , t — ^ 1-t , A — ^ -8 . 

Thus, if x = f ( t ,A,B ) is a solution of (3.6), then 
x = -f(l-t, -B, -A) is also a solution and the 
later solution satisfies the boundary conditions : 

x(O) = -B, x(l) = -A . 

A solution for a given point ( A , B ) generates a solu- 
tion for the reflected point (-B, -A). This reflec- 
ted solution is obtained by the transformation 

x — » -x, t ™-» 1-t 

i.e., a reflection about the t-axis followed by a 
reflection about 't=±' axis. Hence we need 
only consider half of the A,B-plane, i.e., B > -A. 



In the region B < -A, the reflected solution will 
hold. We have the following distinct cases display- 
ing varied nature of the solution. (Kevorkian and 
Cole [58]). The cases are : 

Case - 0. B = A + 1 

Case — 1. A > B— 1 > 0 

Case --2. 0 { J A j < B-l 

Case - 3. B > A+l, ~(B+l) < A < 1-B. 

Case - 4. 3-1 < A < 0, 0<B< A+l 

Case - 5. 0 < B < 1, A > 0. 

Case - 6. B < 0 < A 

First, we cast out the trivial case A = B = 0 ; 
where x = 0 is the solution to (3.6). The two 
outer solutions satisfying the two boundary condi- 
tions respectively are 

x R (t) = t + B - 1 (3.7) 

x R ( t ) = t + A (3.8) 

Case - 0. B = A + 1. 

The two outer solutions are the same, and they 
represent the exact solution. 



Case 


1 . 


A > 3-1 > 0. 


The outer solution x R = t + B - 1 satisfies 
the right boundary condition and the boundary layer 
matched to x R is t 

B l (t) = (B-l) coth (~™"( t+C)) 

T = t/e , 'C = coth" 1 (g^). 

Case - 2. O < |a| < B-l. 

The outer solution x R satisfies the right 
boundary condition and the left boundary layer is 
given by ; 

B l (t) = (B-l) tanh (t+C)] , 

c _ tanh (bIJ) * 

Case -3. B > A+l, -(B+l) < A < 1-B. 

Let t = (1 - A - B)/2. 

The outer solutions x L and x R satisfy both the 
boundary conditions respectively. The inner solu- 
tion at t is : 

I(t}) = tanh (^=~ h) , 

T) = 


(t - t)/e . 


(a shock layer) 



Case ~ 4. 


B~l< A < 0 , 0 < B < A + 1 . 


The outer solution has three pieces i 
x L (t) = t+A 0 < t < -A 

x M (t) = 0 ~A < t < 1-B 

x R (t) = t-l+B 1-B < t < 1 

Here frp provides the match between x r , x,. • 
and f R Q, the match between x,^, x~ • where f R Q 
anc * %vC are ^ wo ex P° nen tial functions satisfying 


lim 
t -»<*> A LC 


f r p("f ) — 0 , 


lim 


— oo L.C 


f (t) = T 


T- f RC< T > - T - T ii“ f RC< T > ■ 0 


Note that f LC and f RC only smoothen the corners 
present in the solution. 


Case - 5. 0 < B < 1, A > 0 

O A £ 

Write f LT = At ;- " 2 s » and f RT » the 

reflection of f L j. In this case, f RT brings the 
solution from A to x ~ 0; f RC brings the solu- 
tion from x = 0 to x R . 

Case — 6 . 

Here, the outer solution is x M = 0$ which is 
joined at the left by f LT and to the right by f RT 
to A > 0 and B < 0 respectively. 



Note that the equation (3.6a) can be written 
in system form as : 



(3.7a) 


ey = x - xy ( 3 . 7 b) 

Ano , then the corresponding linearization is ; 



(3.8a) 
(3.8b) 
(3.8c) 
( 3 . 8d ) 
(3.8e) 


where we assume that we have already initialized 
the functions x Q and y Q . For example, in 
Case 1, we let x (0) = A, x Q (l) = which of- 
course could be taken for all the cases 5 however, 
the numerical experimentation shows that there are 
certain advantages in initializing other functions 
instead of the obvious one, i.e., a straight line 
connecting x Q ( 0 ) = A anci x o^ = 

In the following, we only present the results 
for Case -2, Case -3 and Case-4. Case-0 is omitted 



Mid ~ Minimum of the absolute values of the dif- 
ferences between the asymptotic solution 
and the numerical solution. 


Nit : 


Specifies the iteration number, i.e 
in 'nth iteration'. 



Table 5.1a 


Choice of Stretching Points (Case-2) 


SP : 


e 

Nit 

b l 

bo 


IE -2 

1E~4 

IE «6 

4 

3 

3 

0.036 

0.001 

2E-5 

0.96 

0.996 

99994E-5 


Table 5.1b 


Summary of 

the Comparison 

Results 

(Case-2) 


e 

SS 

C °,b 1 ] 

[b 1 ,b 2 ] 

[b 2 : 

, 1 ] 


NG 

8 

16 

8 


IE -2 

Mad 

24E-3 

88E-3 

28E- 

-4 


Mid 

« — 1 

8 

2E-6 

0E- 

-1 



NG 

16 

16 

8 

IE -4 

Mad 

29E-5 

32E-5 

4E-4 


Mid 

OE-1 

3E-6 

OE-1 


NG 

12 

20 

8 

IE-6 

Mad 

3 IE -6 

15E-6 

IE -5 


Mid 

OE-1 

5E-8 

OE-l 


Table 5.2a 


Choice of Stretching Points (Case-3) 


e 

IE -2 

IE -4 

IE-6 

hit 

3 

4 

5 

SP : b l 

0.46 

0.499 

0.49993 

b 2 

0.522 

0.504 

0.50018 




Table 5.2b 



Summary of 

Comparison Results (Case-3) 


e 

SS 

- 

[o,b 1 ] 

[b x ,b 2 ] 

[b 2 ,l] 


NG 

12 

16 

12 


Mad 

2E-2 

6E-3 

2 IE -3 

IE -2 

Mid 

OE-1 

2E-6 

OE-1 


NG 

12 

16 

12 


Mad 

2E-3 

18E-4 

4E-4 

IE -4 

Mid 

OE-1 

25E-7 

OE-1 


NG 

16 

16 

20 

Mad 

2E-5 

22E-6 

3E-5 

Mid 

OE-1 

6E-8 

OE-1 


JLE-6 




Table 5.3a 


Choice of Stretching Points (Case-4) 



£ 


lt -2 


IE -4 

IE -6 


Nit 

4 


4 

4 


b i 


0.498 


0.4999 

0.49987 

SP : 

b 2 


0.501 


0.5002 

0.50002 


b 3 


0.665 


0.6665 

0.66661 


b 4 


0 . 668 


0.6667 

0.66667 



Table 5.3b 




Summary of 

Comparison 

Results 

» (Case-4) 


e 

ss 

[ 0 ,b x ] 

Cb x ,b 2 ] 

[> 2 

9 ^3 ^ 9 ^4 -I 
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Figure 5 • 1 

Numerical solution to case-2 
-4 

6 =10 , 3 rcJ iteration. 



Figure 5*2 


Numerical solution to case-3 
e=10^ 4 th iteration. 





Figure 5-3 

Numerical solution to case -4 
“4 

€- 10 , 4th iteration. 


5.4 Discussions 


hs almost always it happens, the numerical 
solutions of non-linear problems via iterations is 
a lengthy process. We emphasize on the fact that 
the solution is to be computed after the non-linear 
system has been linearized and not before as e.g. 
in Keller [55], For example, in getting a solution 
for Case-1, e = 10 , we had to solve four linear 

systems obtained in the iteration steps. In spite 
of the long-winding process involved in the numer- 
ical solution, we prefer this approach because of 
the accuracy of the hybrid method suggested in 
Chapter 2 , which we apply for computing the solutions 
of the linear systems. 

When e decreases to zero, the solution approa- 
ches the reduced solution except in the so called 
difficult regions (boundary layers, shock-layers 
and etc.). But the computation of numerical solu- 
tion becomes more difficult when e decreases to zero 
due to the increasing stiffness. However, it may be 
observed from the tables that the numerical solution 
becomes close to the asymptotic solution as e 
decreases. This suggests the robustness of the 
proposed method. Again, neither by analytical nor 



by asymptotical means it is possible to know the 
thickness of the layers. Since, numerically they 
are represented by the stretching sub-intervals, 
this method suggests an approximate thickness of 
the possible layers. But they are by no means 
conclusive. In Case-1, for example, we observe 
that a stretching sub-interval exists around the 
point t = 1, whereas there is no boundary layer 
at this point. This happened because the method 
sensed a possibility of such a layer. In general, 
if one needs finding out whether there is such a 
layer, one should also see whether there is a rapid 
change in the solution in the concerned region. 

From Figure 5.1 it can be seen that there exists no 
such change and hence, no boundary layer near t = 1. 

We note that this method is also applicable to 
nonlinear stiff problems (Kreiss et al. [68] ,Sec.ll ) . 
Stiff problems could also be tackled computationally 
via singular perturbations as e.g. , in Aiken and 
Lapidus [3] and Flaherty and O'Malley [35]. 

We also remark that, once the linearization is 
over, the special methods of Chapter 4 can also be 
applied provided that they are applicable. For 
special non-linearities, such as quasi-linear systems 



of singularly perturbed two-point boundary value 
problems , these methods would probably be appli~ 
cable with advantages (specially the Getting Point 
Technique) . 

Note that the asymptotic solution in Case—4 

has five pieces, where as we observe from Figure 5.3 

that the numerical solution does not show this. 

Instead we have only three pieces. This happens, 

because the asymptotic solution only smoothens' 

1 9 

the corners at t = ^ and ~ *, but, in fact they 
are not necessary, since x^, match at t = ~ 
giving x(l/2) = 0 and also and Xg match 

at t = 2/3 giving x(2/3) = 0. But since there 
are possibilities of corner layers, the grids were 
refined in the neighbourhoods of the points t = 1/2 
2/3. This demonstrates the robustness of the pro- 
posed method. 



Chapter 6 

APPLICATIONS TO OPTIMAL CONTROL 
6.1 Introduction 

The origins of present day control theory can 
be traced back from the feedback theory in 1930's 
by Myguist and Bode and the stability theory of 
nonlinear regulators by Lurie and Krasovski in 1940' s 
The two approaches differ in using frequency domain 
and time domain techniques. By the help of 
Pontryagin's Maximum Principle and Bellman's Dynamic 
Programming the last two and half a decade saw the 
curious merging of these two techniques. Kalman [52] 
synthesized these two methods analytically and 
Wonham [116] expanded the synthesis by the geometric 
methods. Within this period, curious applications 
of singular perturbations to control theory were 
also developed. In fact, long before that control 
engineers were using the simplified models for the 
practical problems they were facing. This new grow- 
ing interest of applying singular perturbation to 
control problems [23] broke them tne happiness that 
what they were doing was legitimate. 



As it happened we did not go to singular per- 
turbation in detail and we do not want to go to 
the singularly perturbed control problems in detail. 
In the last chapter we have suggested and studied 
some numerical methods to solve singularly perturbed 
systems of ordinary differential equations . Our aim 
in this chapter is to apply (some/all of) those 
methods to singularly perturbed optimal control 
problems. Details about these problems and other 
control problems can be had from the review papers 
of Kokotovic et al. [63], O'Malley [83] ,Kokotovic[60] , 
Kokotovic et al. [62] and their references. 

In this chapter we are interested in applying 
the numerical techniques developed in the preceding 
chapters to some optimal control problems. In § 6.2 
the familiar linear state regulators are taken up? 
in & 6.3, the fixed end point optimal control prob- 
lems and in § 6.4, we will consider non-linear 
regulators. § 6.5 gives the observations concerning 
the preceding sections. 

6.2 Linear State Regulators 

The problem here we consider is to select a 
continuous vector function u(t,®) of dimension r 



that minimizes the scalar cost 

J = 2 (z '( t f » e )^(e)z(t f ,e)) 

, t f 

+ 2 S [z 1 (t,e)Q( t,e )z(t,e) 

- 

+ u'(t,e)R(t)u(t,e)]dt ( 2 . 1 a) 

subject to the equations 

x = A 1JL (t,e)x + A 12 (t,e )y + B x (t ,e )u (2.1b) 

ey = A 2J _(t,e )x + A 22 (t,e )y + B 2 (t,s )u (2.1c) 

with the end conditions 

x(t o ,e) = a(e), y(t Q ,e) = P(e ) (2. Id) 

where < t < t < t^ < «>, z = 

x, y are continuously differentiable vector func- 
tions of dimensions m, n respectively? A^, A 12 , 
A 2 i» a 22 , B 1 , Bp, R are smooth matrix functions of 
dimensions mxm, mxn, nxm, nxn, mxr, nxr, rxr res- 
pectively. Also, 

Q u (t,e) Q 12 (t,e) 

Q 22 (t,e) 



Q 








where Q 11 , Q 12 » ^21’ ®22> % i2’ 71 22 are of 

dimensions mxm, mxn, nxm, nxn, mxm, mxn, nxn 

respectively, e > 0 is a small parameter, 
represents d/dt and (') denotes the transpose 
of a matrix; x, y are called, the state variables 
and u, the control. 

In physical problems, the small parameter 
c > 0 often represents negligible time constants, 
small masses, moments of inertia, inductances, 
capacitances etc. For practical purposes, control 
engineers usually neglect the small parameter in 
the model they use for the specific problem and 
then carry on the analysis of the resulting reduced 
model. Of course, they have some reasons to justify 
the approximation! for, in the presence of e , the 
dimensionality of the problem becomes large and also 
the differential equations guiding the states might 
be stiff. But sometimes, it happens that the reduced 
model is a 'bad' approximation to the original model. 
For some specific purposes, it is necessary to carry 
out the analysis without reducing the model (i.e., 



by taking e = 0) . To this end, we employ the 
numerical techniques developed in the preceding 
chapters . 

It is well known that the problem (2.1) has 
a unique solution for fixed s > 0, provided q 
and n are symmetric positive semi definite 
matrices and R is a symmetric positive definite 
matrix. Basically, two different approaches are 
found in the literature that transform (2.1) to 
differential systems. One is via the Hamiltonian 
and the other, via the Riccati transforms. In 
general, for practical purposes, the systems of 
differential equations obtained by either of the 
above approaches, are solved by asymptotic tech- 
niques. We follow the Hamiltonian approach, due 
to its general applicability. 

Introduce the Hamiltonian 

H(x,y ,p,q,u,t, e) = ^ z '^ z + u '^ u ^ + P'^H X + ^ 12 ^ 

+ Bj^u) t q* (A^^x a 22^ ^ 2 ^^ 

( 2 . 2 ) 

for differentiable m,n dimensional vector func- 
tions p(t,s), q(t,<0 (called costate variables) 
that satisfy the equations 



p 


- dH/dx , 


eq 


- dH/dy 


(2.3a) 


with the terminal conditions 

p(t f ,e) - dA/dxjt = t f , eq(t f ,e) = dX/dy] = t f 

(2.3b) 

where a(e) = ~ z’ (t^e )tc(s )z(t^.,e ) , the expres- 
sion which does not come under the integral sign 
in (2.1a). It is known [Pontryagin 1 s theorem , i . e. , 
Ihm 6.4 of Chapter l] that the minimization of 
the cost functional J in (2.1a) 

is equivalent to the minimization of the Hamiltonian 
H in (2.2). Thus dH/du = 0 gives 

u = - R -1 (B^p + B^q). (2.4) 

2 2 

Again, since d H/du = R which is assumed to be 
positive definite, H is minimized by (2.4) and in 
turn J. Hence u is determined by (2.4) provided 
that x,y, p, q are determined as well, x, y, p, q 
are solutions of the following system (q being 
symmetric, Q 2 ^ = Qjp)* 


• 

X 

= A n x 

+ A 12 y 

- B 1 R" 1 B[p - B 1 R“ 1 B^q 

(2.5a) 

© 

8 y 

= A 21 X 

+ A 22 y 

- B 2 R" X B^p - B 2 R" X 3^q 

(2.5b) 

. 

P 

rr X 

- Qj. 2 y 

- - A 21 q 

(2.5c) 



with the end conditions : 


x(t Q ,e ) = a( e) , y(t Q ,e ) = j3(e ) 

(2.5 e ) 

P(tf» £ ) = Tt H x ( tf » £ ) +e ' n 'l2 y ^' t f » e ) 

( 2. 5f ) 

q(t f ,e) = rc[ 2 x(t f ,e) +^ 22 y(t f ,e) 

( 2 . 5g ) 


The two-point boundary-value system (2.5) i s 
a linear homogeneous singularly-perturbed TPBVS ; 
which can be solved by any of the methods suggested 
in the preceding chapters; and. the control u can 
be determined by using equation (2.4). 

Below, we take two familiar examples to illus- 
trate the methods. Solutions to both of them have 
been studied earlier [61,81]. 

The hybrid method is expected to work well due 
to the underlying analysis. So here we will demon- 
strate the other two methods, namely, Boundary Value 
Technique and Cutting Point Technique. 

Example 1 ( [ 8 1 3 ) 

It is required to select the control u(t,e) 
to minimize the cost functional 

j = I[ x 2 (l) + ey 2 (l)j + \ } [ x 2 ( t ) + u 2 (t)]dt (2.6a) 

^ 0 



X 


X 


X = X 

x(0) = 1.3708 

(2.6b) 

ey = -. x - y + u , 

y(0) = 5.9420 

(2.6c) 


Comparing with the equations (2.1) we see that 


1 0 



1 0 

0 • £ 

9 

Q = 

0 0 

> 

b- J 

h J 

II 

h- J 

’ A 12 = 

^21 ^22 ^ 


and &2 — 1* Hence the Hamiltonian reads as follows. 
H( x ,y , p, q , u , t , e) = ™(x^'+u^) + px + q(-x-y+u). 


Minimization of H gives us the following TPBVP : 


x = x 

ey = ~ x 

P = - x 

e q = q 


7 

x(0) 

= 1.3608 

(2.7a) 

y - q 

y(o) 

= 5.9420 

(2.7b) 

p + q , 

p(i) 

= x(l) 

(2.7c) 

9 

q(l) 

= y(i) 

(2.7d) 


0 = u + q 


We want to apply Boundary Value .technique. (2.7a) 
gives the solution of x(t) directly, i.e., 


x ( t ) = 1.3708 exp(t) ( 2 

Eliminating q from (2.7b) - (2.7d), we obtain 



eV = y+(1 '^ (2.8b) 

P + p + 2x - y -ey = 0 (2 . 8c) 

U = x + y + ^Y (2.8d) 

The end conditions are : 

p(l) = 1.3708 e (2.8e) 

ey( 1) = — 1.3708 e (2.8f) 

y( 0 ) = 5.9420. (2.8g) 

00 CO 

Letting p = T p- e 1 , y = 2 y e 1 in (2.8) 

i=o i=o 

and comparing like powers of e, we determine p,y 
upto 0(e) terms, with tj^ = 0.07, 1 2 = 0.96; we 
have (denoting these values with a bar), 

p(t x ) - 2.2065 , p(t 2 ) = 9.5137 (2.9a) 

y(t x ) = 0.2134xl0~ 2 , y(t 2 ) = 1.5327 (2.9b) 

Then we solve the following boundary value problems 

( t = t /e , a = (l-t)/e ) 

dp/dt + e(p - y ~ dy/dt + 2x) = 0 (2.10a) 

d 2 y/dT 2 - e(y + (l-e)x) = 0 


(2.10b) 



y(0) = 5.9420 

(2.10c) 

yUj/e) = 7(t x ), p(t^) = 5(t 1 ) 

(2.10d) 

for 0 < T < t x /e . 


Next , 


•' rip/dcr + e (p ~ y + dy /da + 2x) = 0 

(2.11a) 

d 2 y/da 2 + (y + (l-s)x) = 0 

(2.11b) 

P (0) = 1. 3708e , y(l-t 2 )/e) = y(t 2 ) 

(2.11c) 

( 1+e )y(0) + sdy(0)/da + 1.3708e = 0 

for 0 .< a < (l~t 2 )/s . 

(2. lid) 

Finally , 


• * 

p + p- y~ey + 2x = 0 

(2.12a) 

t 

*< 

i 

t—* 

&) 

X 

n 

o 

(2.12b) 

p(t x ) = p(t 1 ), y(t 1 ) = y(t 1 ) 

(2.12c) 

y(t 2 ) = y(t 2 ) 

(2 . 12d ) 

for < t < t 2 . 



We solve (2.10 - 2.12) with the tolerence limit 
equal to 0.03. 



The exact 
(neglecting the 
given by 


ution of U (t) is approximately 
exponentially small quantities) 


u(t) 


2 1.37 08ft 

3 1+T~ * ex p( ( t~l)/ 8 ) 


(2.13) 


The comparisons of the numerical solution (by 

Boundary Value Technique) and the exact solution 

approximated as in (2.13) are summarized in 
Table 6 . 1 . 


Note that the method worked 
t l» t 2 afresh for different 


even if we do not choose 
e ' s . 


- E -XSffiEi£j2 ([61]) 

Consider a speed control problem for a small 
D/C motor. The state equations are : 


dw/dt = (D/G )i 


eL di/dt = ~ Cw — R i+v 

a 

where w , i and v are speed, current and voltage 
deviations from their nominal values 400 rad/s, 
0.25A and 4.3V. The constants are R = 7.9 q , 

D/G = 1.8654, L = 0.0136H, C = 0.0246 Vs/rad. The 
cost to be minimized is : 

j = + 1 j (w 2 +400i 2 +30v 2 )dt. 



Here % = 





0.01 o 


1 0 

o 0 

, Q = 


J 


o 400 





K 30, v is the control, w and i 
state variables. A u = o, = D/G 

A 22 = " V L ’ B 1 = 0, S 2 = 1/L. 


are the 
A 21 = - C/L, 


Thus, the Hamiltonian is : 


H(w, i, p , q, v, t,e) =I( W 2 + 400i 2 + 30v 2 ) 

+ (D/G)pi + q(- gw i + 1 v ) # 

Minimization of H gives us the following TFBVP : 
(Retaining v instead of q, as = 0 implies 
q = - 30L v); 


dw/rJt = (D/G)i 

(2.14a) 

edi/dt = -(C/L)w _ (R /L)i + (1/L)v 

cl 

(2.14b) 

dp/dt = ~ w ~ 30 Cv 

(2.14c) 

Sdv/dt = (40/3L)i + (D/30LG)p + (R /L)v 

(2.14d) 

with the end conditions : 


w (0) = 400, i( 0) = 0.25 

( 2 . 14 e ) 

p(l) = O.Olw(i), v ( 1 ) = 0. 

(2.14f ) 



We apply Boundary Value Technique to solve (2.14) 
(with tolerence limit = 0.03). The reduced solu- 
tion for the control v is given by 

v(t) = ~ 38.734094 exp(at) + 42.726941 exp(at) 

(2.15) 

with a = 0.49055. 

Comparisons of the numerical solution with the 
reduced solution (2.15) is given in Table 6.1. 

In the table the same terminologies are followed 

Mad, Mid . 

The solution of the linear regulators in 
general, contains boundary layers, whence the left 
boundary layer is stable (in * t * ) and the right 
boundary layer is stable in reverse time (i.e., in 
'-t'). For the occurrence of boundary layers, we 
assume that the system (2.1) is boundary layer 
observable and boundary layer controllable [6l]. 
Further assumptions that we need are : x, y, p, q 
are twice continuously differentiable vector func- 
tions and the conditions needed for the existence 
and uniqueness of the solution to tnis problem. 
Under these assumptions we have obtained a solution 
to the problem using "t — and o— scaled equations. 



The computed results for Example 1 show that 
the solution obtained by the proposed method 
(Boundary Value Technique) approximate the original 
solution well as the small positive parameter 6 
is made smaller, which is a necessary feature for 
numerical methods designed for singular perturba- 
tion problems. From the second example, we observe 
that the solution by the proposed method ( Boundary 
Value Technique) does approach the reduced solution 
as e approaches zero. (In the second example, we 
have taken the reduced solution for comparison 
purposes, for the stiffness involved in the 
problem) . 


It has been observed that the method works 
well on linear optimal control problems. The 
applicability of the present methods to fixed-end. 
point problems and non-linear regulators come as 
possible generalizations. In later sections, these 
possibilities will be investigated. 



Table 6.1 


Summary of the Comparison Results for the Examples 



Ex 

- 1 


Ex - 

2 

s 

| £ ' ' 


iwad 

0.711E-3 


Mad 

0.099561 

X JU *** 1 

i 'iid 

0.35E-6 

0.5 

Mid 

0.028142 


Mad 

0.225E-3 


Mad 

0.10271 

IE -2 

Mid 

0.1E-6 

0.1 

Mid 

0.000928 


Mad 

0. 155E-3 


Mad 

0.013471 

IE-3 

Mid 

O.OE-1 

0.05 

Mid 

0.000021 


Mad 

0.478E-9 


Mad 

0.000173 

IE -4 

Mid 

O.OE-1 

0.01 

Mid 

0.000003 


Mad 

0.248E-11 


Mad 

0.000002 

IE -5 



0.001 




Mid 

O.OE-1 


Mid 

0.000000 




6.3 Singularly Perturbed Linear Fixed End-Point 
Optimal Control and Cheap Control 

In § 6.2 we have considered the linear 
regulators and demonstrated the efficiency of the 
numerical method (Boundary value technique) for 
such problems. A natural extension is the fixed 
end-point problems, i.e,, when the end-points 
x(t Q ), x(t f ), y(t 0 ), y(t f ) of the (state) trajec- 
tories are fixed, (see equation - (2. Id)). .Pre- 
cisely, the problem is to select control u(t,e) 
that minimizes the cost (integral cost) 


J = 

, tf 

— / ( z ’ z + u'Ru)dt 

2 b 

(3.1a) 

subject to the state 


• 

X = 

A U X + A 12 y + B lY 

(3.1b) 

ey = 

A 21 x + A 22 y + B 2 y 

(3.1c) 

with 

the boundary conditions 


x(t 0 , 

e) = x°, y(t 0 ,e) = Y° 

(3. Id) 

X ( t j: f 

e) = x f , y(t f ,e) = y f 

(3.1e) 


where z = + C 2 (t)y , g, C 2 being 

given matrices of functions of t; and other things 



reading the same as defined in § 6.2. 

Applying Pontryagin's Principle (Theorem 6.4 
of Chapter 1) to (3.1) we find that the optimal 
trajectory, (x,y) the co-states (p,q) and the 
control u are given as 


u 

= - R~ 1 (B|p + 

B^q) 

(3.2a) 

■» 

X 

= A n x + A 12 y 

- S n p - S 12 q 

(3.2b) 

• 

ey 

= m 21 x + A 22 y 

" S i2 p " S 22 q 

(3.2c) 

* 

P 

= qc x x ~ c* 

C ? y - A{ lP - 

(3. 2d) 

eq 

= - C^C 1 x - C' 2 

C 2 y - A’ jP - A^q 

(3.2e) 


with boundary conditions (3. Id - 3.1e) and 

S . . = B . R" 1 B , i,j = 1,2. 

1J i 3 

For the applications of this problem one is 
referred to Sannuti [102], Khalil [59], Hadlock [40] 
It is known that this problem is not robust to 
singular perturbations. In fact, the left-boundary 
layer is stable in time whereas the right boundary 
layer is stable in reverse-time (negative time) [37, 
40, 102]. This is usually solved either by a Riccat 
equations approach or by the method of Matched Asym- 

The difficulties associated with 


ptotic expansions. 



these approaches are with finding positive and 
negative definite solutions of the Riccati equa- 
tions [114] or with 'matching' in a hypothetical 
oomain. The numerical methods developed in the 
previous chapters apply well to this case. The 
IP6VP contained in (3,2) with the boundary condi- 
tions (3. Id ~ 3.1e) can be solved in a similar 
manner as the corresponding problem in $ 6.2. 


Now we consider another related linear control 
problem, namely Cheap Control. Consider a cost 
functional (to be minimized) 

^f 

J(e) =1 / (x'Qx + e^u'Ru)dt (3.3a) 

' y 

U, being symmetric positive semidefinite and R, 
symmetric positive definite smooth matrix functions 
respectively. The control constraints being in the 
time-*varying system form 


x = Ax + Su 


t < t < t f 
0 X 


(3.3b) 


with an initial condition 

,, \ o (3.3 

x(t Q ) = x 

where A(t), B(t) are smooth matrix functions of 
dimensions nxn and nxr ; x and u are column 



vectors of dimensions n and r respectively? 

£ 0 is a small parameter. (Note that the appear™ 

ance of e is in the cost and not in the states). 
Classical control theory implies that [52] the 
optimal control 

u ••= ~ £ 2 R -1 B'p (3.4a) 

where the co-state variable p satisfies 
p = - Qx ~ A'p (3.4b) 

and the end-condition 

p(t f ) = 0. (3.4c) 

Eliminating u, from (3.3 - 3.4), we have the 
following singularly perturbed TP3VP : 


s 2 Ax -• BR" 

' X B ' p 

(3.5a) 

-Qx - A’p 


(3.5b) 

*°, 

p(t f ) = o. 

(3.5c) 


which can be solved numerically. However, as in 
O'Malley and Jameson [87] we will consider a special 
subclass of such cheap control problems. 


e-L' by requiring that 



BIQB. = 0 


j = 0 , 1 ,.. . ,L-2 


H ! 

L- 


0 


on 


with 


b j 

= - 

- 

j 2 i ? 


The special subclass refers to the problems that 
belong to Case-L, for some L. 

It is shown in [87] that for Case-L, the boundary 
layer variables ( t and c previously) must be 

chosen by 



1/L 

a = (t f -t)/(e) (3.6) 


Hence, when applying the special methods (i.e., 
Boundary value technique and cutting point technique), 
the boundary layers must be stretched according to 

( 3 . 6 ). 

Note that we need not transform the problems into 
an n+Lr - dimensional problem as in [87] for having 
a numerical solution. This suggests an obvious 
advantage of the suggested numerical approaches over 
those employed in the afore-cited reference. For 
more details, one is advised to see appropriate 
references in Kokotovic [60], 



6.4 Nonlinear State Regulators 

In this section we consider an optimization 
problem of Lagrange type where the system equations 
axe stiff and the system can be expressed in terms 
of a singular perturbation parameter. Specifically, 
we consider the problem of minimizing the scalar 
cost 

/ *f 

J - n(x(t f ), y(t f ), e) + / cp(x,y,u,t,e)dt (4.1a) 

t o 

subject to the state constraints 

x = f ( x ,y , u,t ,e ) (4.1b) 

ey = g(x,y,u,t,e ) (4.1c) 

x(t Q ) = x°, y(t Q ) = y° (4. Id) 

where x(t,e) and y(t,e) are smooth state vari- 
ables of dimensions m and n respectively, 
u(t,e) is the control vector of dimension r, the 
real number s is a small positive parameter and 
' denotes d/dt. It is desired to select the 
control u to minimize J and thereby determine 
the state trajectory. 

In physical problems, the small positive 
parameter s represents inductances, capacitances. 



moments of inertia, and small masses, which control 
engineers often neglect. Instead of taking the 
problem (4.1), they consider the reduced problem 
(i.e. by letting e= 0) because of the difficul- 
ties that arise from the high dimensionality and 
the stiffness of the system (4.1). But in many 
cases , the reduced model becomes a 'bad' approx- 
imation to the original model and it loses many 
of the important properties. 

When s is close to zero, the problems are 
such that near t = t and t = there exists 

0 i 

no uniformly valid solution. In other words, the 
solution of the original problem does not converge 
uniformly to the solution of the reduced one as 
£ — > 0. Such a situation is commonly referred to 
as the occurrence of a boundary layer at t = t Q 
and at t = t^. Mostly, the solutions to these 
problems have been constructed with the help of 
matched asymptotic expansions (cf. Kokotovic [60]). 
Then the solutions of the boundary-layer equations 
and that of the reduced equations are combined in 
some way to give an approximate solution over the 
entire interval. As usual, define the Hamiltonian 
H with the help of co-state variables p and q 



of dimensions m and n respectively : 

H(x, y , P > q > u , t , e ) = cp + p'f + q * g 

here ( ) represents transpose of a vector or 
matrix. The necessary conditions for optimality 
follow from Pontryagin's principles [Thm. 6.4 of 
Chapter l]. In addition to (4.1b) - (4. Id), we 

have 

p = -AH/dx , eq = -dH/dy (4.2a) 

p(t f ,e) = dw(t f ,e)/dx, eq(t-,s) = drc(t-,e)/dy 

1 f (4.2b) 

dli/ hu = 0 , (4.2c) 

2 2 

d H/0u“ > 0 is sufficient for the minimization 
of J . Now, our task is to solve (4.1b) - (4.2b) 
numerically and determine u from (4.2c). Note 
that if we want to use the Hybrid method we have 
to follow the linearization suggested in Chapter 5. 

On the other hand, iterations can be avoided by 
using the special methods developed in Chapter 4. 

In the following, Cutting point Technique is modified 
to solve the above problem class. Boundary Value 
Technique can also be modified similarly. However, 
if the problem has a turning point, one is advised 
to follow the Hybrid method strictly. 

For the purpose of solving the TPBVP (4.1b-4.2b), 



we assume the following : 

(A~l) The strong Legendre-Clebsch condition 
d 2 H/du 2 >0 holds. 

(A-2) tt(x,y, e ) = e(x, e y, e ) for some functions 

p (., .). 

(A-3) The reduced problem of (4.1b) ~ (4.2b) has 

a unique solution (with appropriate boundary- 
conditions) 

( A-4 ) The boundary layers are locally stable 

(with dichotomies, i.e., including reverse- 
stability near t = t^) and there is no 
turning point except in the neighbourhoods 
of t Q and t f (stretching factor = e ). 

Assumption A-l guarantees that H u = 0 in (4.2c) 
can be solved uniquely to yield u = r)(x,y,p,q,t ,e ) 
(see Thm. 6.5 of Chapter l); A-2 tells that ti is a 
slowly varying function of the fast state variable 
y; A -3 will help in finding out values of x,y,p,q 
at the points t x and t 2 (necessary for the Cut- 
ting Point Technique). Assumption A-4 forms the 
basis of scaling the equations by T and cr. With 
these assumptions (of which A-3 and A-4 are exclu- 
sively necessary for the special method : i.e., 
Cutting point Technique), we proceed to describe 
the technique in a stepwise manner, as follows : 



S tep 1 


Choose ^ close to t Q and t^ 

respectively . We expect that the boundary layers 
lie within the intervals [t^tj and [t 2 ,t f ], (as 
experience suggests we let the thickness of the 
boundary layers be 0< | ©In s | ) ) . With this expecta- 
tion, find the values of x,y,p,q at t 1 and. t 2 
from the reduced equations corresponding to (4.1b - 
4.2b), neglecting the boundary conditions for y(t ) 
and q(t^)‘, i.e., we evaluate x(t^), y(t^) etc. 
from the following either by analytical or by 
numerical means : 

x ~ f, p = ~dH/dx, g = 0 = dH/dy (4.3a) 

x(t Q ) = x°, p(t f ) = d7i(t f ,0)/dx (4.3b) 

1 2 

Denote x(t^) = x » x(t 2 ) = x , etc. 


S tep , 2 . 


In this step, our airn is to get an approximate 
solution to the problem (4.1b) - (4.2b). To this eno 

we lot 


(t - t n )/e 


(ti - tJA 


<3 == (t f - t)/ & (4.4a) 
a 2 = (t f - t 2 )/e (4.4b) 


and form the r-scaled 


and the a-scaled equation of 



(4.1b) -- (4.2b), using the transformations (4.4a). 
be supply the t - scaled, original and cr-scaled 
equations with the boundary conditions as in the 

following . • ■ ■ 

t — scaled Original o— scaled 

x,y at t Q x,y at t-^ x,y at 

p,q at t^ p,q at tj p,q at t^ . 

Let X -- <x,y,p,q> denote the stacked (2m+2n) -vector 
( x ' ,y ' ,p' , q * ) * and (note the use of < ., ., ., .,>) 


I 

0 


o 

0 , 

'm+n 

m+n 

B = 

m+n 

m+n 

0 

m+n 

°m+n 

y ! 

i 

1 

. 

^m+n 

"'‘m+n 


where is the identity matrix of order k and 

0 . is the square zero matrix of order k, and pro- 

K 

ceed as follows : 


Solve (4.5) for t £ [0, \] : 

X = < f, g, ~ dH/dx, -dH/dy > 

AX(O) + BX(\) =<x°,y°,p 1 ,q 1 > 
Solve (4.6) for t £ [t J _,t 2 ] J 
X = < f, e'- l g, -dH/dx, -e 1 dH/dy > 

ax ( t x ) + ex(t 2 ) = < x\y\p 2 ,<i 2 > 


(4.5a) 

(4.5b) 


(4.6a) 


(4.6b) 



(4.7a) 


Finally, solve (4.7) for a 6 [O,^] J 

* 

^ ^ f -g , sdH/dx, dli/dy > 

‘~' / ‘ ; (°) + a '^* 2 ) =<x2 »y 2 . d7t(t f ,e)/dx, 

e" i d7i(t f ,e )/dy > (4.7b) 

we combine the solutions of (4.5 - 4.7) 
using the inverse transformations of (4.4a) : 

t = ex + t 0 on [0, xj and 

t == t f ~ sc on [0 ,a 2 ], 

btep 3 

In this step, we describe briefly how to solve 
the nonlinear TPBVP's obtained above. Here we apply 
centered difference scheme. Our problems are of 

the following form : 

NX - X - F(t,X) = 0, AX( a) + BX(b) = C (4.8) 

whore N is a nonlinear operator and C is a 
2(m+n ) -dimensional constant vector. To ensure the 
convergence of the Newton's iteration scheme (des- 
cribed below) we assume that F satisfies the 
following conditions : 

(t,X) is continuous for, i»j = 1,2, . . ,2m+2n . 

(i) ^ 



How , with some initial guess Z Q (we might take 
the reduced solution) we compute the sequence 
^ s> > s = 0,1,.., by the following : 

S.+l = Z s + Z s (4.12a) 

where Zg is the solution of the linear algebraic 

system 

•|f (^s ) "s + *(2 S ) = 0 (4.12b) 

Note that the error in this procedure is 0 (h). 

(cf. [55, § 3.3]). 

Step 4 

In this step, our aim is to confirm whether 
our choice of t^ and t 2 is acceptable, or 
whether we need to seek a better choice. To this 
end, we fix up a tolerence limit on the difference 
between the derivatives (given below) computed from 
the three different regions [t 0 ,t 1 ], [t 1 ,t 2 J, [t 2 ,t f ] 
at the points t^ t,,. So, we compute X(tp^ from 
(4.5; _ 4.6) and denote them by X_(r 1 ) and ^ + (tp 



respectively; we also compute X(t 2 ) from (4.6 - 

4.7) to denote them by Xjy and X + (t 9 ) 
respectively. We then suggest the following cri- 
terion which is to be satisfied by our computed 

solution ,* 

I U + ( tj ) ~ < 6 , j = 1,2 (4.13) 

where r, > 0 is some prescribed tolerence limit. 

If (4.13) is satisfied, we stop, else, we choose 
different t^ greater than,and t 2 less than that 
previously chosen and repeat the steps 1-4. 

in the following, we take two examples to 
illustrate the method. The first example is a 
quasi linear one taken from Sannuti [103] and the 
second is from Freedman and Granoff [36] (also 
Ardema [4]) which is a nonlinear. 


Examp l e 3 

Consider the problem of minimizing 


J = A / ( x 2 + u 2 )dt 
z o 

with the state equations 


y + e x 


x(0) 


o 

X 


11.45 


ey 


- x - y + u, 


y(0) = y° = 1.985 • 



The Hamiltonian is 


! ! ( x , y , p , q , u , t , e ) = - ■|(x 2 +u 2 ) + p(y+ £ x 2 ) + q(-x-y+u). 


The ITjjVP corresponding 

to (4.1b) - (4 

,2b) is : 

* 2 
x ~ y + e x , 

X 

II 

X 

o 

(4.14a) 

m 

c y - ~ x - y + q , 

y(0) = y° 

(4.14b) 

p ~ - x ■ 2exp + q, 

o 

II 

* — I 

"cu 

(4.14c) 

* 

q ■= - p + q , 

o 

II 

> — { 

cr 

(4.14d) 

u ~ q . 


(4.15) 

The reduced solution of 

(4.14) without 

the condi- 

tions y(O) = y° and q(l) = 0 is ; 


x = x°(l - t/2) , 

y = - x°/2 

, (4.16a) 

p = x c (l ~ t)/2 , 

q = x°(l - 

t)/2.(4.16b) 

'fe choose t^ and t 2 

corresponding to various 

e's as follows : 



e 0.1 

0.001 

0.00001 

t l 0,05 

0.005 

0.00005 

t 2 0.05 

0.006 

0.00C05 

Next we compute x(tj) 

etc. from (4.16) with 



those choices of t i - i o . 

3 1,2 , and denote these by 

“ sJ ’ 8 E « x >y>P.<l 1 . J = 1,2. The problems 
in step 2 to be solved are the following : 


Deno t 

e t 

= t/e , 

c = (1- 

~t)/e 


= t x / £ 

; and 

n 2 " 

U 

-tp)/ e , 







The ^ 

equations 

corresponding 

to (4 

.5) are : 

dx/di 

- 

2 2 

£ X + e y 


9 

x(0) = 

= x° 

(4.17a) 

dy/di 

as 

- x - y + 

q 

9 

y (o) = 

= y° 

(4.17b) 

H p/dt 

=: 

- e x - 2e 

2 xp + q 

9 

p( T 1 ) 

= P X 

(4.17c) 

Hq/ rjf 

::*S 

- p + q 


9 

q( T, ) 

= q 1 

(4.17d 


Similarly the equations corresponding to (4.6 - 4.7) 

might be formed. Note that all r-scaled, original 
and o-scaled equations corresponding to (4.5) - (4.7) 
satisfy four assumptions in step 3. Finally, we 
combine the results of (4.17). &c. ) to get a solu- 

tion for (4.14), using the transformations t = ex 
or t = 1 ~ so in the appropriate domains. Since 
in this modification we have followed the conventional 
approach of taking finite difference schemes first 
and apply linearization afterwards, we give the full 
numerical results in the following tables instead of 



jui,t the comparison results as we were following 
I r * 'V i o' is Ly , The numerical solutions are given in 

TnMes (6.2) - (6.4). 

hxam r. le 4 

Select the controls u and v to minimize 


the terminal cost 
J = x( l,e ) + - ey(l 

cons traints 

x = y + u - tv + e 

. 2 
ey - - y - tu + 2v + 

x(O f e ) = x° = 28.239, 

Considering the corre: 

n - J ) 

H “ pf + qg , 

we have the following 

x = f , ey = g, 

2„ 

p= - ep-eQ, 
e q = _p+q”2eyq> 


e) subject to the state 

x+u) = f 
s 2 (x+y 2 ) = g 

y(0,e) = y° = 5.984. 
ponding Hamiltonian (with 

IPBVP : 

x(0) = x° , y(O) = y° 

p(l, e) = 1 

/ \ 1 
qU*«' = 2 ’ 

v = tp/4q, 


u = tq/2(p~e) , 



whose reduced solution is : 


3 *3 3 9 

■' = x o “ 8 1 * y = - 8 1 » p = q = 1 • 

V. choose t^ and t ,2 for corresponding e’s as 

follows : 


c 


O.l 

0.001 

0.00001 

n 


u.05 

0.005 

0.0005 

el 


0.05 

0.005 

0.0005 

and compute w (tj ) = 

= w"" 1 and solve the following 

systems 

* 

(see equations (4.5) - 

- (4.7)) for 

t e [ 

0 # 

V : 



rix/ifl 

= 

, 2 
e (y + u 

- tv) + e 2 (x H 

h u) 

dy/di 

- 

-(y + tu) 

9 2 

+ 2 v + e (x 

+ y 2 ) 

dp/ dt 

SS 

- e 2 (p + < 

sq) 


dq/ di 

sr 

- p + q " 

2 e 2 yq 


x ( 0 ) = 

X° 

, y(o) = 

y° » p( t i) = 

p 1 , q(x x ) = q 1 

f or t 

e 

[ti,t 2 ] : 




y + u 2 - tv + e (x + u) 

,(y + tu - 2v 2 ) + s 2 ( x + y ) 


p = - . e ( p + e g ) 
eq = -p + q- 2 e 2 yq . 

X(t x ) = X 1 , y(t x ) = V 1 - p 1 ^ 



q(t 2 ) = cf 2 • 



f iTiril I.y , for o e [o, 0 2 ] : 

C >£(y + U 2 - tV ) - £ 2 ( X + u ) 

'*'/ / f!, “ ~ y + tu - 2v 2 - e?(x + y 2 ) 

rift/'": <T “ e 2 (p + eq) 

Hq/-;C = P - q + 2 £?yq 

x(°2 } s y^ 2 } = Y 2 ’ P(°) = 1, q(0) = 1/2. 

Mote that all the three TPBVP's above satisfy the 

conditions in step 3. 

I text , we combine the solutions of the above 
TPiiVi 1 g with the help of the transformations 
t - ct and t = 1 - ec in the corresponding domains. 
The numerical solutions of this example are given in 
Tables (6.5) - (6.7) along with the asymptotic solu- 
tions for different values of e . 

In both of the examples we have taken the tole- 
rence limit T = 0.0005 whence the requirement (4.13) 
is satisfied immediately. However, for specific 
problems one might need an experienced eye to mini- 
mize wasting machine-time in choosing appropriate 
t } and t 2 . In the tables the results for the most 
sensitive state y(t,0 out of the trajectories (x,y) 


is given 



Table 6.2 


numerical Results for Example 3 ( e = o.l) 


0.0 

0 . 0( 500c 1 
0 . O0UG05 
O.UOOOI 
0.00005 
0 . UUQ5 
c.oOl 
o.oub 
0 . 1*1 
0.05 
0.1 

0 . 3 
0.4 

o.;> 

0.6 

0.7 

0.8 

0.9 

0.95 

1.0 


h = 1/50 

1.985000 

2.125632 

2.228346 

2.302382 

2.932187 

3.201035 

3.312389 

3.021392 

2.8510869 

- 0.972318 

- 3.922875 

- 6.337021 

- 6.638210 

- 6.189023 

- 5.547271 

- 4.893127 

- 4.288932 

- 3.721013 

- 2.923105 

- 2.221032 

0.000829 


Y (t,s) 
h = 1/100 

1.985000 

2.210835 

2.293701 

2.808126 

3.123460 

3.570831 

3.705723 

3.454023 

2.851453 

- 0.975937 

- 3.922258 

- 6.336842 

- 6.631289 

- 6.189000 

- 5.545323 

- 5.894350 

- 4.287721 

- 3.697003 

- 2.998929 

- 1.902587 

0.000829 


Asymptotic 

solution 

4.098675 

4.098545 

4.098024 

4.097373 

4.092168 

4.033754 

3.959164 

3.454186 

3.851472 

- 0.975958 

- 3.922261 

- 6.336836 

- 6.631239 

- 6.189005 

- 5.545291 

- 5.894299 

- 4.287563 

- 3.697209 

- 2.998204 

- 2.515795 

- 1.845201 


Table 6.3 

I'-'miQTital Results for Example 3 ( e = 0.001) 


t 



1 h = 1/50 

( ; . i 

1.985000 

{ . 1. w'U-U < • 1 

1.723109 

0 . Cv. : OCC.J5 

1.689022 

L . i K. u. ] 

1.653125 

i. , ( 't '5 

1.392370 

Iv. €05 

-3.169983 

i) ,i)Ol 

-6.300897 

( . » : 6 

*-11.013457 

(, ■ * u 1. 

-11.020893 


.*10.396324 

C. 1 

-9.662935 

0.2 

-8.333954 

0. 3 

-7.171729 

0.4 

-6.157028 

0.5 

-5.260932 

u . 6 

-4.477029 

0.7 

-3.781232 

0.8 

-3.160902 

0.9 

-2.595230 

0.9b 

-2.000210 

o 

* 

<HI 

0,000229 


Y (t,c) 


h = 1/100 

Asymptotic 

solution 

1.985000 

2.007126 

1.823192 

1.993972 

1.801236 

1.941489 

1.792793 

1.876179 

1.369999 

1.365296 

-3.170833 

-3.170497 

-6.309312 

-6.309256 

-11.012239 

-11.012114 

-11.020203 

-11.020199 

-10.395291 

-10.395398 

-9.660323 

-9.660406 

-8.332809 

-8.332531 

-7.172530 

-7.172223 

-6.156093 

-6.156116 

-5.263920 

-5.263758 

-4.477192 

-4.477197 

-3.780632 

-3.780617 

-3.160202 

-3.160013 

-2.600001 

-2.602911 

-1.078923 

-2.344632 

0.000222 

-2.084680 


Table 6.4 


cal Results for Example 3 ( e = O.OOOOl) 

| y ( t . s) 

! h - 1/50 h = 1/100 Asymptotic 

i solution 


>G< n ; 1 


1 . 905000 


O.i * l 005 
U, i 'trt t, .1 
0 .-05 

O. t >01 
S . t o‘> 

O, ul 
u. u5 
0. 1 
0.2 
0.3 

t ) • *1 

0.5 
0 . 6 

o,7 

0.8 

0.9 

0.95 

1.0 


1.000213 
- 2.990703 
- 5.002317 
- 11.085792 
- 11.166782 
- 11.155299 
- 11.092005 
- 11.011927 
- 10.387729 
- 9.653991 
- 8.326953 
- 7.157729 
- 6.150009 
- 5,255952 
- 4.402937 
- 3.707621 
- 3.000293 
- 2.008232 
- 1.780236 
0.000002 


1.985000 

0.932124 

-3.002037 

-5.302938 

-11,083092 

-11.164023 

-11.156228 

-11.091487 

-11.011072 

-10.386792 

-9.653029 

-8.326871 

-7.167592 

-6.150999 

-5.259032 

-4.203921 

-3.774029 

-2.892316 

-1.999997 

-1.002138 

- 1.000000 


1.986211 

0.734018 

-3.191253 

-3.331528 

-11.083111 

-11.164464 

-11.156343 

-11.091569 

-11.011099 

-10.383838 

-9.653043 

-8.326885 

-7.167548 

-6.151807 

-5.259311 

-4.472181 

-3.774646 

-3.152732 

-2.593980 

-2.334749 

-2.087062 



Table 6.5 


Numerical Results for Example 4 (s = 0.1) 

y(t,e) 


1/100 


Asymptotic 

solution 


0.0 

0.000001 
0. 000005 
0.00001 
0.00005 
0.0005 
0.001 
0. 005 
0.01 
0.05 
0.1 
0.2 
0.3 
0.4 
0.5 
0.6 
0.7 
0.8 
0.9 
0.95 
1.0 


h = 1/50 
5.984000 
5.899231 
5.723450 
5.529322 
5.203124 
5.598346 
5.802341 
5.690235 
5.497831 
3.823092 
2.272346 
0.832937 
0.326821 
0.124592 
0.033023 
-0.018235 
-0.062029 
-0.092135 
-0.113215 
-0.094921 
-0.009009 


5.984000 

5.821023 

5.692319 

5.402132 

5.011295 

5.492358 

5.892031 

5.709235 

5.448921 

3.689234 

2.236002 

0.844792 

0.326602 

0.124492 

0.033921 

-0.019623 

-0.062788 

-0.098899 

-0.114499 

-0.095507 

0.000209 


5.008991 
5.008931 
5.008692 
5.008394 
5 . 006006 
5.979208 
5.949574 
5.717763 
5.440750 
3.659778 
2.235121 
0.844804 
0.326596 
0. 124476 
0.033984 
-0.019716 
-0.062713 
-0.099565 
-0.114521 
-0.095578 
-0.021886 



Table 6.6 

Numerical Results for Example 4 ( e = O.OOl) 


t 

y(t,e) 

h = 1/50 

h = 1/100 

Asymptotic 

solution 

o 

d 

5 . 984000 

5.984000 

5.984250 

0.000001 

5.978321 

5.978293 

5.978269 

0 . 000005 

5.958220 

5.954889 

5.954405 

0.00001 

5.924002 

5.924799 

5.924708 

0.00005 

5.693836 

5.692025 

5.692407 

0.0005 

3.633275 

3.629008 

3.629730 

0.001 

2.212301 

2.202012 

2.201641 

0.005 

0.041927 

0.040678 

0.040567 

O 

• 

o 

H ' 

0.000998 

0.000402 

0 . C 00497 

0.05 

- 0.000987 

- 0.000682 

- 0.000625 

0.1 

- 0.003828 

- 0.003397 

- 0.003375 

0.2 

- 0.018235 

- 0.014689 

- 0.014500 

00 

• 

o 

- 0.032729 

- 0.033295 

- 0.033125 

0.4 

- 0.059927 

- 0.059292 

- 0.059250 

0.5 

- 0.095273 

- 0.092792 

- 0.092875 

0.6 

- 0.129985 

- 0.134995 

- 0.134000 

0.7 

- 0.183829 

- 0.182692 

- 0.182625 

0.8 

- 0.230956 

- 0.238012 

- 0.238750 

0.9 

- 0.300125 

- 0.302257 

- 0.302375 

0.95 

- 0.123792 

- 0.321785 

- 0.337000 

1.0 

- 0.000095 

- 0.120937 

- 0.123972 





Table 6.7 


Numerical Results for Example 4 (e = O.OOOOl) 


t 

y(t,e) 

h = 1/50 

h = 1/100 

Asymptotic 

solution 

0.0 

5.984000 

5.984000 

5.984003 

0.000001 

5.415832 

5.414228 

5.414550 

0.00000b 

3.628821 

3.629823 

3.629482 

0.00001 

2.202586 

2.201872 

3.201393 

0.00005 

0.048201 

0.040923 

0.040322 

0.0005 

0.002357 

0.000089 

0.000002 

0.001 

0.000092 

0.000000 

0. 000002 

0.005 

-0.000001 

-0.000004 

-0.000007 

0.01 

-0.000098 

-0.000032 

-0.000035 

0.05 

-0.002015 

-0.000996 

-0.000934 

0.1 

-0.003892 

-0.003788 

-0.003746 

0.2 

-0.014026 

-0.014823 

-0.014995 

0.3 

-0.033082 

-0.033789 

-0. 033744 

0.4 

-0.059087 

-0.059987 

-0.059993 

0.5 

-0.093000 

-0.093788 

-0.093741 

0.6 

-0.138820 

-0.134028 

-0.134990 

0.7 

-0.182992 

-0.183599 

-0. 183739 

0.8 

-0.238892 

-0.239888 

-0.239988 

0.9 

-0.303557 

-0.303786 

-0.303736 

0.95 

-0.338119 

-0.338400 

-0. 338423 

1.0 

-0.099267 

-0.100278 

-0. 124989 



6.5 Discussions 


In § 6.2 •- 6.4, we have indicated how the 
methods developed in the previous chapters might be 
used lor solving optimal control problems numeri- 
cally. We have concentrated on the special methods 
suggested in Chapter 4. The Hybrid method is 
expected to behave well in this application area. 

The main feature of the Hamiltonian systems 
that arise in connection with optimal control prob- 
lems is that they contain two boundary layers, in 
general, at the two ends of the interval of inte- 
gration. The left boundary layer is stable in 
time whereas the right boundary layer is stable 
in negative time, the thickness being of O(jelnej) 
or roughly 0(e)? (see the assumptions in $ 4.2) 
which we exploit in choosing the stretchings 
t = (t - t )/e and a = (t f - t)/e at the left 
and right boundary layers. For the linear problems 
we assume that the systems are boundary layer obs- 
ervable and boundary layer controlable which guarentee 
the existence of the boundary layers and in turn, 
justifies the stretchings. Under these crucial 
assumptions (see [60-61]), we have obtained solu- 
tions to the linear optimal control problems using 
x and o-scaled equations. 



The computed results for examples for linear 
regulators ( § 6.2) show that the solution obtained 
by the proposed method approximates the original 
solution better as the mesh is made finer and the 
parameter e is made smaller. In the second 
example, we have taken the reduced solution for 
comparison purposes, because of the stiffness in- 
volved in the problem. Examples are not given for 
the fixed end point problems and cheap control 
problems ( § 6.3) for the obvious similarity of 

these problems with the linear regulators. For 
non-linear regulators however examples are given 
to illustrate the methods. We have taken the 
rapidly varying trajectory component y(t,e) for 
comparison purposes. None of the examples has 
right boundary-layer behaviour. When we decrease 
e to zero, and when the 

mesh size h is made finer, the convergence of the 
numerical solution- to the reduced one becomes 
faster. We have given the numerical solutions 
obtained by the proposed methods along with those by 
the asymptotic method, in the tables ( § 6.4) 
hoping that it would facilitate comparison. 

We also find that the solutions are close 
enough except at the boundaries, where the solution 



obtained by the proposed method is obviously better 
than the other. The same dichotomies regarding 
the stability of the boundary layers still hold 
for non-linear problems. Our assumptions guaran- 
tee that there is no turning point in the middle, 
whence we are justified in trying to find a solu- 
tion in the boundary layers with the stretchings t 
and a. Since the thickness of the boundary layers 
is not yet determined with analytical or asymptotic 
means for the Hamiltonian systems, and asymptotic 
methods available to solve this type of problem 
almost avoid this question, it is suggested that 
we try to get a solution by the proposed methods 
where we would be able to know more accurately 
about the thickness of the boundary layers (rather 
than just 0(e) or O(jeinej). These methods have 
also the advantage of simplicity over the non-tri- 
vial matching procedures for such problems. 



Chapter 7 


CONCLUSION 


The wide applicability of singular perturba- 
tion on the one hand and the inherent challenges 
of singular perturbation on the other stimulated 
an untiring interest in the mathematicians. It 
is the latter that attracts every person with a 
mathematical bent. However, the reason for nume- 
rical analysts to poke their noses involves both. 

A numerical analyst can never be satisfied to know 
only the nature of a solution of a particular 
class of problem(s). He is also interested in get- 
ting a solution concretely (quantitatively ) .A pri- 
mary interest of him is an efficient method. He 
requires both conceptual and computational effi- 
ciency. The method should, according to a nume- 
rical analyst, be simple, efficient and accurate. 

The first and foremost criterion is accuracy. To 
be a accurate, a numerical method must be convergent 
(i.e., non-accumulation of errors and etc.) and the 
numerical solution that it produces must possess 
all the qualitative properties. 



In trying to find out numerical methods for 
singularly perturbed two-point boundary-value 
problems, we have implicitly narrowed down our 
interest. This limitation of scopes of the 
applicability of the methods was guided by two 
essential properties : the boundary-layer pheno- 
mena and oscillation. We have considered problems 
whose solutions show boundary layer behaviours by 
excluding the possibility that the eigen-values 
of the system matrix be purely imaginary and 
eradicating oscillations by smoothness conditions. 
Singularly perturbed systems having boundary layer 
behaviours always possess two modes : slow state 
and fast state (in control theory language.). The 
fast state again consists of two modes : rapidly 
growing and rapidly decaying. These properties 
are projected well in decoupling the slow and fast 
modes in Chapter 2. The decoupling transformation 
introduced has the advantage of simplicity over the 
usual Lypunov-type transformation The analytic 
decoupling (compared to the decoupling of differ- 
ence schemes) facilitates applying Euler's 
schemes to the fast mode. However, a similar prob- 
lem remains there such as decoupling the difference 
equations (by an implicit multi-step method) by 



using Lypunov-type transformations [62] (parallel 
to iviatheij [74]). We emphasize that an implicit 
scheme rather than an explicit scheme is to be 
used, due to the failure of Lypunov-type transfor- 
mations on the difference equations yielded by the 
explicit one-step method, (see § 2.4 - 2.5). 

Guided by simplicity and efficiency, we have 
also suggested some special methods to handle some 
sub-classes of the singular perturbation problems, 
namely, the problems without turning points. Though 
we have only seen those special methods to work 
well for homogeneous problems, they are expected 
to work well also for non~homogeneous ones. For 
non-linear problems, the special methods have the 
advantage that they can be applied with the usual 
process of linearization. In contrast the general 
purpose method, we have developed in Chapter 2 needs 
the linearization of a non-linear problem before 
applying the difference schemes. This was achieved 
by approximating the non-linear function F(t,z), 

(in z = F (t,z)) by a sequence of linear functions 
(see Chapter 5). 

The numerical experimentations carried out for 
the reliability and realizability of the methods 



developed report on the affirmative. It was obs- 
erved that the methods approximate the exact/ 
asymptotic solutions with less labour. When the 
parameter e approaches zero, the numerical sol- 
ution also approaches the reduced solution (pre- 
serving the non-uniformity characteristic of such 
problems). It was also observed that the expected 
accuracy could always be achieved without much 
effort. 

There are at least two ways to combat stiff- 
ness involved in the computation of solutions of 
singularly perturbed systems. One is to design 
a better computer (hardware) and the other is to 
design a better algorithm (soft ware). These two 
are not necessarily mutually exclusive. However, 
our emphasis was on the latter alternative. ... 

All the computations were carried out on 
DEC-1090 computer system available at Indian 
Institute of Technology, Kanpur. 
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